cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
prt_header.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*PrtHeader print program's header, including luminosities and ionization parameters */
4 #include "cddefines.h"
5 #include "phycon.h"
6 #include "iso.h"
7 #include "rfield.h"
8 #include "radius.h"
9 #include "called.h"
10 #include "thermal.h"
11 #include "dense.h"
12 #include "continuum.h"
13 #include "ipoint.h"
14 #include "prt.h"
15 #include "freebound.h"
16 
17 void PrtHeader(void)
18 {
19  long int i,
20  ip2500,
21  ip2kev;
22  double absbol,
23  absv,
24  alfox,
25  avpow,
26  bolc,
27  gpowl,
28  pballog,
29  pionl,
30  qballog,
31  qgaml,
32  qheiil,
33  qhel,
34  ql,
35  qxl,
36  radpwl,
37  ratio,
38  solar,
39  tcomp,
40  tradio;
41 
42  DEBUG_ENTRY( "PrtHeader()" );
43 
44  if( !called.lgTalk )
45  {
46  return;
47  }
48 
49  if( prt.lgPrtCitations )
50  {
51  /* the print citations command */
53  fprintf( ioQQQ, "\n" );
55  fprintf( ioQQQ, "\n\n" );
56  }
57 
58  fprintf( ioQQQ, " %4ldCellPeak",rfield.nflux );
60  fprintf( ioQQQ, " Lo");
61  fprintf( ioQQQ,PrintEfmt("%9.2e", rfield.anumin(0) ));
62  fprintf( ioQQQ, "=%7.2fm Hi-Con:", 1e-10*RYDLAM/rfield.anumin(0) );
64  fprintf(ioQQQ," Ryd E(hi):");
66  fprintf(ioQQQ,"Ryd E(hi): %9.2f MeV\n", rfield.egamry()*1e-6*EVRYD );
67 
68  if( prt.xpow > 0. )
69  {
70  prt.xpow = (realnum)(log10(prt.xpow) + radius.pirsq);
71  qxl = log10(prt.qx) + radius.pirsq;
72  }
73  else
74  {
75  prt.xpow = 0.;
76  qxl = 0.;
77  }
78 
79  if( prt.powion > 0. )
80  {
81  pionl = log10(prt.powion) + radius.pirsq;
82  avpow = prt.powion/rfield.qhtot/EN1RYD;
83  }
84  else
85  {
86  pionl = 0.;
87  avpow = 0.;
88  }
89 
90  /* >>chng 97 mar 18, break these two out here, so that returns zero
91  * when no radiation - had been done in the print statement so pirsq was printed */
92  if( prt.pbal > 0. )
93  {
94  pballog = log10(MAX2(1e-30,prt.pbal)) + radius.pirsq;
95  qballog = log10(MAX2(1e-30,rfield.qbal)) + radius.pirsq;
96  }
97  else
98  {
99  pballog = 0.;
100  qballog = 0.;
101  }
102 
103  if( radius.pirsq > 0. )
104  {
105  fprintf( ioQQQ, " L(nu>1ryd):%9.4f Average nu:",pionl);
106  PrintE93(ioQQQ, avpow);
107  fprintf( ioQQQ," L( X-ray):%9.4f L(BalC):%9.4f Q(Balmer C):%9.4f\n",
108  prt.xpow, pballog, qballog );
109  if( pionl > 47. )
110  {
111  fprintf(ioQQQ,"\n\n WARNING - the continuum has a luminosity %.2e times greater than the sun.\n",
112  exp10( pionl-log10(SOLAR_LUMINOSITY) ) );
113  fprintf(ioQQQ," WARNING - Is this correct? Check the luminosity commands.\n\n\n");
114  }
115  }
116  else
117  {
118  fprintf( ioQQQ, " I(nu>1ryd):%9.4f Average nu:",pionl);
119  PrintE93(ioQQQ, avpow);
120  fprintf( ioQQQ," I( X-ray):%9.4f I(BalC):%9.4f Phi(BalmrC):%9.4f\n",
121  prt.xpow,
122  log10(MAX2(1e-30,prt.pbal)),
123  log10(MAX2(1e-30,rfield.qbal)) );
124  }
125 
126  if( rfield.qhe > 0. )
127  {
128  qhel = log10(rfield.qhe) + radius.pirsq;
129  }
130  else
131  {
132  qhel = 0.;
133  }
134 
135  if( rfield.qheii > 0. )
136  {
137  qheiil = log10(rfield.qheii) + radius.pirsq;
138  }
139  else
140  {
141  qheiil = 0.;
142  }
143 
144  if( prt.q > 0. )
145  {
146  ql = log10(prt.q) + radius.pirsq;
147  }
148  else
149  {
150  ql = 0.;
151  }
152 
153  if( radius.pirsq != 0. )
154  {
155  fprintf( ioQQQ,
156  " Q(1.0-1.8):%9.4f Q(1.8-4.0):%9.4f Q(4.0-20):"
157  "%9.4f Q(20--):%9.4f Ion pht flx:",
158  ql,
159  qhel,
160  qheiil,
161  qxl);
163  }
164  else
165  {
166  fprintf( ioQQQ,
167  " phi(1.0-1.8):%7.4f phi(1.8-4.0):%7.3f phi(4.0-20):"
168  "%7.3f phi(20--):%7.3f Ion pht flx:",
169  ql,
170  qhel,
171  qheiil,
172  qxl );
174  }
175  fprintf( ioQQQ,"\n");
176 
177  /* estimate alpha (o-x) */
178  if( rfield.anu(rfield.nflux-1) > 150. )
179  {
180  /* the ratio of fluxes is nearly 403.3^alpha ox */
181  ip2kev = ipoint(147.);
182  ip2500 = ipoint(0.3645);
183  if( rfield.flux[0][ip2500-1] > 1e-28 )
184  {
185  ratio = (rfield.flux[0][ip2kev-1]*rfield.anu(ip2kev-1)/rfield.widflx(ip2kev-1))/
186  (rfield.flux[0][ip2500-1]*rfield.anu(ip2500-1)/rfield.widflx(ip2500-1));
187  }
188  else
189  {
190  ratio = 0.;
191  }
192 
193  if( ratio > 0. )
194  {
195  alfox = log(ratio)/log(rfield.anu(ip2kev-1)/rfield.anu(ip2500-1));
196  }
197  else
198  {
199  alfox = 0.;
200  }
201  }
202  else
203  {
204  alfox = 0.;
205  }
206 
207  if( prt.GammaLumin > 0. )
208  {
209  gpowl = log10(prt.GammaLumin) + radius.pirsq;
210  qgaml = log10(prt.qgam) + radius.pirsq;
211  }
212  else
213  {
214  gpowl = 0.;
215  qgaml = 0.;
216  }
217 
218  if( prt.pradio > 0. )
219  {
220  radpwl = log10(prt.pradio) + radius.pirsq;
221  }
222  else
223  {
224  radpwl = 0.;
225  }
226 
227  if( radius.pirsq > 0. )
228  {
229  fprintf( ioQQQ,
230  " L(gam ray):%9.4f Q(gam ray):%9.4f L(Infred):%9.4f Alf(ox):%9.4f Total lumin:%9.4f\n",
231  gpowl, qgaml, radpwl, alfox, log10(MAX2(SMALLFLOAT,continuum.TotalLumin)) +
232  radius.pirsq );
233  }
234  else
235  {
236  fprintf( ioQQQ,
237  " I(gam ray):%9.4f phi(gam r):%9.4f I(Infred):%9.4f Alf(ox):%9.4f Total inten:%9.4f\n",
238  gpowl, qgaml, radpwl, alfox, log10(MAX2(SMALLFLOAT,continuum.TotalLumin)) +
239  radius.pirsq );
240  }
241 
242  /* magnitudes */
243  if( radius.lgPredLumin )
244  {
245  solar = absbol = -38.;
247  {
248  solar = log10(continuum.TotalLumin) + radius.pirsq - log10(SOLAR_LUMINOSITY);
249  absbol = 2.5*(log10(MBOL_ZERO_LUMINOSITY/SOLAR_LUMINOSITY) - solar);
250  }
251 
252  /* absv = 4.79 - 2.5 * (LOG10(MAX(1e-30,continuum.fluxv)) + pirsq - 18.742 -
253  * 1 0.016)
254  * allen 76, page 197 */
255  absv = -38.;
257  absv = -2.5*(log10(MAX2(1e-30,continuum.fluxv)) + radius.pirsq - 20.654202);
258 
259  /* >>chng 97 mar 18, following branch so zero returned when no radiation at all */
260  if( continuum.fbeta > 0. )
261  {
262  continuum.fbeta = (realnum)(log10(MAX2(1e-37,continuum.fbeta)) + radius.pirsq);
263  }
264  else
265  {
266  continuum.fbeta = 0.;
267  }
268 
269  bolc = absbol - absv;
270  fprintf( ioQQQ,
271  " log L/Lsun:%9.4f Abs bol mg:%9.4f Abs V mag:%9.4f Bol cor:%9.4f nuFnu(Bbet):%9.4f\n",
272  solar,
273  absbol,
274  absv,
275  bolc,
276  continuum.fbeta );
277  }
278 
279  rfield.cmcool = 0.;
280  rfield.cmheat = 0.;
281 
282  for( i=0; i < rfield.nflux; i++ )
283  {
284  /* CSIGC is Tarter expression times ANU(I)*3.858E-25 */
286  rfield.csigc[i];
287 
288  /* Compton heating with correction for induced Compton scattering
289  * CSIGH is Tarter expression times ANU(I)**2 * 3.858E-25 */
291  rfield.csigh[i]*(1. + rfield.OccNumbIncidCont[i]);
292  }
293 
294  /* all of following need factor of 10^-15 to be true rates */
295  rfield.cmheat *= dense.eden*1e-15;
296  rfield.cmcool *= dense.eden*4./TE1RYD*1e-15;
297 
298  if( rfield.cmcool > 0. )
299  {
300  rfield.lgComUndr = false;
301  tcomp = rfield.cmheat/rfield.cmcool;
302  }
303  else
304  {
305  /* underflow if Compt cooling rate is zero */
306  rfield.lgComUndr = true;
307  tcomp = 0.;
308  }
309 
310  /* check whether energy density temp is greater than compton temp */
311  if( phycon.TEnerDen > (1.05*tcomp) )
312  {
313  thermal.lgEdnGTcm = true;
314  }
315  else
316  {
317  thermal.lgEdnGTcm = false;
318  }
319 
320  /* say some ionization parameters */
321  fprintf( ioQQQ, " U(1.0----):");
322  PrintE93( ioQQQ, rfield.uh);
323  fprintf( ioQQQ, " U(4.0----):");
325  fprintf( ioQQQ, " T(En-Den):");
327  fprintf( ioQQQ, " T(Comp):");
328  PrintE93( ioQQQ,tcomp );
329  fprintf( ioQQQ, " nuJnu(912A):");
331  fprintf( ioQQQ, "\n");
332 
333  /* some occupation numbers */
334  fprintf( ioQQQ, " Occ(FarIR):");
336  fprintf( ioQQQ, " Occ(H n=6):");
337 
338  if( iso_sp[ipH_LIKE][ipHYDROGEN].n_HighestResolved_local + iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local >= 6 )
339  {
340  long ip6p = iso_sp[ipH_LIKE][ipHYDROGEN].QuantumNumbers2Index[6][1][2];
341  PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ip6p].ipIsoLevNIonCon-1]);
342  }
343  else
344  {
345  PrintE93( ioQQQ, 0.);
346  }
347  fprintf( ioQQQ, " Occ(1Ryd):");
349  fprintf( ioQQQ, " Occ(4R):");
351  fprintf( ioQQQ, " Occ (Nu-hi):");
353  fprintf( ioQQQ, "\n");
354 
355  /* now print brightness temps */
356  tradio = rfield.OccNumbIncidCont[0]*TE1RYD*rfield.anu(0);
357  fprintf( ioQQQ, " Tbr(FarIR):");
358  PrintE93( ioQQQ, tradio);
359  fprintf( ioQQQ, " Tbr(H n=6):");
360  if( iso_sp[ipH_LIKE][ipHYDROGEN].n_HighestResolved_local + iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local >= 6 )
361  {
362  long ip6p = iso_sp[ipH_LIKE][ipHYDROGEN].QuantumNumbers2Index[6][1][2];
363  PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ip6p].ipIsoLevNIonCon-1]*TE1RYD*rfield.anu(iso_sp[ipH_LIKE][ipHYDROGEN].fb[ip6p].ipIsoLevNIonCon-1));
364  }
365  else
366  {
367  PrintE93( ioQQQ, 0.);
368  }
369  fprintf( ioQQQ, " Tbr(1Ryd):");
370  PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1]*TE1RYD*rfield.anu(iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1));
371  fprintf( ioQQQ, " Tbr(4R):");
372  PrintE93( ioQQQ, rfield.OccNumbIncidCont[iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1]*TE1RYD*rfield.anu(iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].ipIsoLevNIonCon-1));
373  fprintf( ioQQQ, " Tbr (Nu-hi):");
375  fprintf( ioQQQ, "\n");
376 
377  if( tradio > 1e9 )
378  {
379  fprintf( ioQQQ,
380  " >>>The radio brightness temperature is very large, %.2eK at %.2ecm. Is this physical???\n",
381  tradio, 1e-8*RYDLAM/rfield.anu(0) );
382  }
383 
384  /* skip a line */
385  fprintf( ioQQQ, " \n" );
386  return;
387 }
double TEnerDen
Definition: phycon.h:108
realnum * csigh
Definition: rfield.h:271
realnum q
Definition: prt.h:280
t_thermal thermal
Definition: thermal.cpp:6
void PrintE93(FILE *, double)
Definition: service.cpp:1019
double exp10(double x)
Definition: cddefines.h:1383
double widflx(size_t i) const
Definition: mesh.h:147
realnum qbal
Definition: rfield.h:339
realnum ** flux
Definition: rfield.h:70
const realnum SMALLFLOAT
Definition: cpu.h:246
realnum * outlin_noplot
Definition: rfield.h:191
#define PrintEfmt(F, V)
Definition: cddefines.h:1364
realnum qhe
Definition: rfield.h:339
t_phycon phycon
Definition: phycon.cpp:6
realnum fluxv
Definition: continuum.h:80
realnum ** outlin
Definition: rfield.h:191
realnum pradio
Definition: prt.h:280
FILE * ioQQQ
Definition: cddefines.cpp:7
bool lgTalk
Definition: called.h:12
bool lgPrtCitations
Definition: prt.h:294
vector< freeBound > fb
Definition: iso.h:481
void CloudyPrintReference()
Definition: service.cpp:1890
double anu(size_t i) const
Definition: mesh.h:111
t_dense dense
Definition: global.cpp:15
realnum fx1ryd
Definition: prt.h:280
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:15
realnum pbal
Definition: prt.h:280
const int ipH1s
Definition: iso.h:29
t_continuum continuum
Definition: continuum.cpp:6
realnum qhtot
Definition: rfield.h:339
double cmcool
Definition: rfield.h:277
t_rfield rfield
Definition: rfield.cpp:9
realnum * ConInterOut
Definition: rfield.h:156
float realnum
Definition: cddefines.h:124
realnum uh
Definition: rfield.h:347
realnum GammaLumin
Definition: prt.h:289
realnum uheii
Definition: rfield.h:350
double cmheat
Definition: rfield.h:277
void PrtHeader(void)
Definition: prt_header.cpp:17
realnum * OccNumbIncidCont
Definition: rfield.h:119
t_radius radius
Definition: radius.cpp:5
multi_arr< long, 3 > QuantumNumbers2Index
Definition: iso.h:490
double anumin(size_t i) const
Definition: mesh.h:139
t_prt prt
Definition: prt.cpp:14
realnum pirsq
Definition: radius.h:148
bool lgPredLumin
Definition: radius.h:144
double TotalLumin
Definition: continuum.h:71
realnum xpow
Definition: prt.h:280
realnum fbeta
Definition: continuum.h:81
const int ipH_LIKE
Definition: iso.h:64
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
const int ipHELIUM
Definition: cddefines.h:349
realnum qgam
Definition: prt.h:280
realnum powion
Definition: prt.h:280
realnum qheii
Definition: rfield.h:339
realnum * csigc
Definition: rfield.h:271
double egamry() const
Definition: mesh.h:88
void DatabasePrintReference()
Definition: service.cpp:1907
double eden
Definition: dense.h:201
#define MAX2(a, b)
Definition: cddefines.h:828
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
double anumax(size_t i) const
Definition: mesh.h:143
bool lgEdnGTcm
Definition: thermal.h:77
void PrintE82(FILE *, double)
Definition: service.cpp:920
realnum qx
Definition: prt.h:280
const int ipHYDROGEN
Definition: cddefines.h:348
bool lgComUndr
Definition: rfield.h:284
long int nflux
Definition: rfield.h:48
t_called called
Definition: called.cpp:4
long int ipeak
Definition: prt.h:288