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