cloudy  trunk
 All Data Structures 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-2008 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 "physconst.h"
6 #include "phycon.h"
7 #include "iso.h"
8 #include "rfield.h"
9 #include "radius.h"
10 #include "called.h"
11 #include "thermal.h"
12 #include "dense.h"
13 #include "continuum.h"
14 #include "ipoint.h"
15 #include "prt.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  fprintf( ioQQQ, " %4ldCellPeak",rfield.nflux );
51  fprintf( ioQQQ, " Lo");
52  fprintf( ioQQQ,PrintEfmt("%9.2e", rfield.anu[0] - rfield.widflx[0]/2. ));
53  fprintf( ioQQQ, "=%6.2fcm Hi-Con:",
54  9.117e-6/(rfield.anu[0] - rfield.widflx[0]/2.) );
56  fprintf(ioQQQ," Ryd E(hi):");
58  fprintf(ioQQQ,"Ryd E(hi): %9.2f MeV\n", rfield.egamry*0.0000136 );
59 
60  if( prt.xpow > 0. )
61  {
62  prt.xpow = (realnum)(log10(prt.xpow) + radius.pirsq);
63  qxl = log10(prt.qx) + radius.pirsq;
64  }
65  else
66  {
67  prt.xpow = 0.;
68  qxl = 0.;
69  }
70 
71  if( prt.powion > 0. )
72  {
73  pionl = log10(prt.powion) + radius.pirsq;
74  avpow = prt.powion/rfield.qhtot/EN1RYD;
75  }
76  else
77  {
78  pionl = 0.;
79  avpow = 0.;
80  }
81 
82  /* >>chng 97 mar 18, break these two out here, so that returns zero
83  * when no radiation - had been done in the print statement so pirsq was printed */
84  if( prt.pbal > 0. )
85  {
86  pballog = log10(MAX2(1e-30,prt.pbal)) + radius.pirsq;
87  qballog = log10(MAX2(1e-30,rfield.qbal)) + radius.pirsq;
88  }
89  else
90  {
91  pballog = 0.;
92  qballog = 0.;
93  }
94 
95  if( radius.pirsq > 0. )
96  {
97  fprintf( ioQQQ, " L(nu>1ryd):%9.4f Average nu:",pionl);
98  PrintE93(ioQQQ, avpow);
99  fprintf( ioQQQ," L( X-ray):%9.4f L(BalC):%9.4f Q(Balmer C):%9.4f\n",
100  prt.xpow, pballog, qballog );
101  if( pionl > 47. )
102  {
103  fprintf(ioQQQ,"\n\n WARNING - the continuum has a luminosity %.2e times greater than the sun.\n",
104  pow( 10. , pionl-log10(SOLAR_LUMINOSITY) ) );
105  fprintf(ioQQQ," WARNING - Is this correct? Check the luminosity commands.\n\n\n");
106  }
107  }
108  else
109  {
110  fprintf( ioQQQ, " I(nu>1ryd):%9.4f Average nu:",pionl);
111  PrintE93(ioQQQ, avpow);
112  fprintf( ioQQQ," I( X-ray):%9.4f I(BalC):%9.4f Phi(BalmrC):%9.4f\n",
113  prt.xpow,
114  log10(MAX2(1e-30,prt.pbal)),
115  log10(MAX2(1e-30,rfield.qbal)) );
116  }
117 
118  if( rfield.qhe > 0. )
119  {
120  qhel = log10(rfield.qhe) + radius.pirsq;
121  }
122  else
123  {
124  qhel = 0.;
125  }
126 
127  if( rfield.qheii > 0. )
128  {
129  qheiil = log10(rfield.qheii) + radius.pirsq;
130  }
131  else
132  {
133  qheiil = 0.;
134  }
135 
136  if( prt.q > 0. )
137  {
138  ql = log10(prt.q) + radius.pirsq;
139  }
140  else
141  {
142  ql = 0.;
143  }
144 
145  if( radius.pirsq != 0. )
146  {
147  fprintf( ioQQQ,
148  " Q(1.0-1.8):%9.4f Q(1.8-4.0):%9.4f Q(4.0-20):"
149  "%9.4f Q(20--):%9.4f Ion pht flx:",
150  ql,
151  qhel,
152  qheiil,
153  qxl);
155  }
156  else
157  {
158  fprintf( ioQQQ,
159  " phi(1.0-1.8):%7.4f phi(1.8-4.0):%7.3f phi(4.0-20):"
160  "%7.3f phi(20--):%7.3f Ion pht flx:",
161  ql,
162  qhel,
163  qheiil,
164  qxl );
166  }
167  fprintf( ioQQQ,"\n");
168 
169  /* estimate alpha (o-x) */
170  if( rfield.anu[rfield.nflux-1] > 150. )
171  {
172  /* the ratio of fluxes is nearly 403.3^alpha ox */
173  ip2kev = ipoint(147.);
174  ip2500 = ipoint(0.3645);
175  if( rfield.flux[ip2500-1] > 1e-28 )
176  {
177  ratio = (rfield.flux[ip2kev-1]*rfield.anu[ip2kev-1]/rfield.widflx[ip2kev-1])/
178  (rfield.flux[ip2500-1]*rfield.anu[ip2500-1]/rfield.widflx[ip2500-1]);
179  }
180  else
181  {
182  ratio = 0.;
183  }
184 
185  if( ratio > 0. )
186  {
187  alfox = log(ratio)/log(rfield.anu[ip2kev-1]/rfield.anu[ip2500-1]);
188  }
189  else
190  {
191  alfox = 0.;
192  }
193  }
194  else
195  {
196  alfox = 0.;
197  }
198 
199  if( prt.GammaLumin > 0. )
200  {
201  gpowl = log10(prt.GammaLumin) + radius.pirsq;
202  qgaml = log10(prt.qgam) + radius.pirsq;
203  }
204  else
205  {
206  gpowl = 0.;
207  qgaml = 0.;
208  }
209 
210  if( prt.pradio > 0. )
211  {
212  radpwl = log10(prt.pradio) + radius.pirsq;
213  }
214  else
215  {
216  radpwl = 0.;
217  }
218 
219  if( radius.pirsq > 0. )
220  {
221  fprintf( ioQQQ,
222  " L(gam ray):%9.4f Q(gam ray):%9.4f L(Infred):%9.4f Alf(ox):%9.4f Total lumin:%9.4f\n",
223  gpowl, qgaml, radpwl, alfox, log10(continuum.TotalLumin) +
224  radius.pirsq );
225  }
226  else
227  {
228  fprintf( ioQQQ,
229  " I(gam ray):%9.4f phi(gam r):%9.4f I(Infred):%9.4f Alf(ox):%9.4f Total inten:%9.4f\n",
230  gpowl, qgaml, radpwl, alfox, log10(continuum.TotalLumin) +
231  radius.pirsq );
232  }
233 
234  /* magnitudes */
235  if( radius.lgRadiusKnown )
236  {
237  solar = log10(continuum.TotalLumin) + radius.pirsq - 33.5828;
238  absbol = 4.75 - 2.5*solar;
239 
240  /* absv = 4.79 - 2.5 * (LOG10(MAX(1e-30,continuum.fluxv)) + pirsq - 18.742 -
241  * 1 0.016)
242  * allen 76, page 197 */
243  absv = -2.5*(log10(MAX2(1e-30,continuum.fluxv)) + radius.pirsq - 20.654202);
244 
245  /* >>chng 97 mar 18, following branch so zero returned when no radiation at all */
246  if( continuum.fbeta > 0. )
247  {
248  continuum.fbeta = (realnum)(log10(MAX2(1e-37,continuum.fbeta)) + radius.pirsq);
249  }
250  else
251  {
252  continuum.fbeta = 0.;
253  }
254 
255  bolc = absbol - absv;
256  fprintf( ioQQQ,
257  " log L/Lsun:%9.4f Abs bol mg:%9.4f Abs V mag:%9.4f Bol cor:%9.4f nuFnu(Bbet):%9.4f\n",
258  solar,
259  absbol,
260  absv,
261  bolc,
262  continuum.fbeta );
263  }
264 
265  rfield.cmcool = 0.;
266  rfield.cmheat = 0.;
267 
268  for( i=0; i < rfield.nflux; i++ )
269  {
270  /* CSIGC is Tarter expression times ANU(I)*3.858E-25 */
272  rfield.csigc[i];
273 
274  /* Compton heating with correction for induced Compton scattering
275  * CSIGH is Tarter expression times ANU(I)**2 * 3.858E-25 */
277  rfield.csigh[i]*(1. + rfield.OccNumbIncidCont[i]);
278  }
279 
280  rfield.cmheat *= dense.eden*1e-15;
281 
282  /* 6.338E-6 is k in inf mass Rydbergs */
283  rfield.cmcool *= dense.eden*4.*6.338e-6*1e-15;
284 
285  /* all of following need factor of 10^-15 to be true rates */
286  if( rfield.cmcool > 0. )
287  {
288  rfield.lgComUndr = false;
289  tcomp = rfield.cmheat/rfield.cmcool;
290  }
291  else
292  {
293  /* underflow if Compt cooling rate is zero */
294  rfield.lgComUndr = true;
295  tcomp = 0.;
296  }
297 
298  /* check whether energy density temp is greater than compton temp */
299  if( phycon.TEnerDen > (1.05*tcomp) )
300  {
301  thermal.lgEdnGTcm = true;
302  }
303  else
304  {
305  thermal.lgEdnGTcm = false;
306  }
307 
308  /* say some ionization parameters */
309  fprintf( ioQQQ, " U(1.0----):");
310  PrintE93( ioQQQ, rfield.uh);
311  fprintf( ioQQQ, " U(4.0----):");
313  fprintf( ioQQQ, " T(En-Den):");
315  fprintf( ioQQQ, " T(Comp):");
316  PrintE93( ioQQQ,tcomp );
317  fprintf( ioQQQ, " nuJnu(912A):");
319  fprintf( ioQQQ, "\n");
320 
321  /* some occupation numbers */
322  fprintf( ioQQQ, " Occ(FarIR):");
324  fprintf( ioQQQ, " Occ(H n=6):");
325 
327  {
328  long ipH6p = iso.QuantumNumbers2Index[ipH_LIKE][ipHYDROGEN][6][1][2];
330  }
331  else
332  {
333  PrintE93( ioQQQ, 0.);
334  }
335  fprintf( ioQQQ, " Occ(1Ryd):");
337  fprintf( ioQQQ, " Occ(4R):");
339  fprintf( ioQQQ, " Occ (Nu-hi):");
341  fprintf( ioQQQ, "\n");
342 
343  /* now print brightness temps */
344  tradio = rfield.OccNumbIncidCont[0]*TE1RYD*rfield.anu[0];
345  fprintf( ioQQQ, " Tbr(FarIR):");
346  PrintE93( ioQQQ, tradio);
347  fprintf( ioQQQ, " Tbr(H n=6):");
349  {
350  long ipH6p = iso.QuantumNumbers2Index[ipH_LIKE][ipHYDROGEN][6][1][2];
352  }
353  else
354  {
355  PrintE93( ioQQQ, 0.);
356  }
357  fprintf( ioQQQ, " Tbr(1Ryd):");
359  fprintf( ioQQQ, " Tbr(4R):");
361  fprintf( ioQQQ, " Tbr (Nu-hi):");
363  fprintf( ioQQQ, "\n");
364 
365  if( tradio > 1e9 )
366  {
367  fprintf( ioQQQ,
368  " >>>The radio brightness temperature is very large,%10.2eK at%10.2ecm. Is this physical???\n",
369  tradio, 9.115e-6/rfield.anu[0] );
370  }
371 
372  /* skip a line */
373  fprintf( ioQQQ, " \n" );
374  return;
375 }

Generated for cloudy by doxygen 1.8.3.1