cloudy trunk
Loading...
Searching...
No Matches
prt_zone.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/*PrtZone print out individual zone results, call by iter_end_check at very
4 * end of zone calculations */
5#include "cddefines.h"
6#include "physconst.h"
7#include "iso.h"
8#include "grainvar.h"
9#include "pressure.h"
10#include "wind.h"
11#include "conv.h"
12#include "trace.h"
13#include "magnetic.h"
14#include "dense.h"
15#include "called.h"
16#include "dynamics.h"
17#include "h2.h"
18#include "mole.h"
19#include "secondaries.h"
20#include "opacity.h"
21#include "colden.h"
22#include "geometry.h"
23#include "hmi.h"
24#include "rfield.h"
25#include "thermal.h"
26#include "radius.h"
27#include "phycon.h"
28#include "abund.h"
29#include "hydrogenic.h"
30#include "ionbal.h"
31#include "elementnames.h"
32#include "atomfeii.h"
33#include "prt.h"
34#include "taulines.h"
35
36void PrtZone(void)
37{
38 char chField7[8];
39 char chLet,
40 chQHMark;
41 long int i,
42 ishift,
43 nelem ,
44 mol;
45 double hatmic;
46
47 DEBUG_ENTRY( "PrtZone()" );
48
49 if( thermal.lgUnstable )
50 {
51 chLet = 'u';
52 }
53 else
54 {
55 chLet = ' ';
56 }
57
58 /* middle of zone for printing
59 rmidle = radius.Radius - radius.drad*0.5*radius.dRadSign;
60 dmidle = radius.depth - radius.drad*0.5; */
61
62 /* option to print single line when quiet but tracing convergence
63 * with "trace convergence" command */
64 if( called.lgTalk || trace.nTrConvg )
65 {
66 /* print either ####123 or ###1234 */
67 if( nzone <= 999 )
68 {
69 sprintf( chField7, "####%3ld", nzone );
70 }
71 else
72 {
73 sprintf( chField7, "###%4ld", nzone );
74 }
75
76 fprintf(ioQQQ, " %7.7s %cTe:",chField7, chLet);
78 fprintf(ioQQQ," Hden:");
79 PrintE93(ioQQQ,dense.gas_phase[ipHYDROGEN]);
80 fprintf(ioQQQ," Ne:");
81 PrintE93(ioQQQ,dense.eden);
82 fprintf(ioQQQ," R:");
83 PrintE93(ioQQQ,radius.Radius_mid_zone );
84 fprintf(ioQQQ," R-R0:");
85 PrintE93(ioQQQ,radius.depth_mid_zone);
86 fprintf(ioQQQ," dR:");
87 PrintE93(ioQQQ,radius.drad);
88 fprintf(ioQQQ," NTR:%3ld Htot:",conv.nPres2Ioniz);
90 fprintf(ioQQQ," T912:");
91 fprintf(ioQQQ,PrintEfmt("%9.2e",opac.TauAbsGeo[0][iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] ));
92 fprintf(ioQQQ,"###\n");
93
94 if( trace.nTrConvg )
95 {
96 fprintf( ioQQQ, " H:%.2e %.2e 2H2/H: %.2e He: %.2e %.2e %.2e\n",
97 dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
98 dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
99 2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN],
100 dense.xIonDense[ipHELIUM][0]/SDIV(dense.gas_phase[ipHELIUM]),
101 dense.xIonDense[ipHELIUM][1]/SDIV(dense.gas_phase[ipHELIUM]),
102 dense.xIonDense[ipHELIUM][2]/SDIV(dense.gas_phase[ipHELIUM])
103 );
104 }
105 }
106
107 /* now return if not talking */
108 if( !called.lgTalk && !trace.nTrConvg )
109 {
110 return;
111 }
112
113 /* lgDenFlucOn set to true in zero, only false when variable abundances are on,
114 * lgAbTaON set true when element table used */
115 if( !dense.lgDenFlucOn || abund.lgAbTaON )
116 {
117 fprintf( ioQQQ, " Abun:" );
118 for( i=0; i < LIMELM; i++ )
119 {
120 fprintf( ioQQQ,PrintEfmt("%8.1e", dense.gas_phase[i] ));
121 }
122 fprintf( ioQQQ, "\n" );
123 }
124
125 /*-------------------------------------------------
126 * print wind parameters if windy model */
127 if( !wind.lgStatic() )
128 {
129 double fac;
130 /* find denominator for fractional contributions */
131 if( wind.AccelTotalOutward == 0. )
132 fac = 1.;
133 else
134 fac = wind.AccelTotalOutward;
135 fprintf( ioQQQ,
136 " Dynamics wind V:%.3e km/s a(grav):%.2e a(tot):%.2e Fr(cont):%6.3f "
137 "Fr(line):%6.3f \n",
138 wind.windv/1e5 ,
139 -wind.AccelGravity,
140 wind.AccelTotalOutward,
141 wind.AccelCont/ fac,
142 wind.AccelLine/fac );
143
144 /* print advection information */
145 if( dynamics.lgAdvection )
146 DynaPrtZone();
147 }
148
149 /* print line with radiation pressure if significant */
150 if( pressure.pbeta > .05 )
152
153 /*---------------------------------------------------- */
154
155 hatmic = 0.;
157 for(mol = 0; mol < mole_global.num_calc; mol++) {
158 if (mole_global.list[mol]->parentLabel.empty() && mole_global.list[mol]->nAtom.find(elHydrogen) !=
159 mole_global.list[mol]->nAtom.end())
160 hatmic += mole.species[mol].den*mole_global.list[mol]->nAtom[elHydrogen];
161 }
162 ASSERT(hatmic > 0.);
163 hatmic = (dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHYDROGEN][1])/hatmic;
164
165 fprintf( ioQQQ, " Hydrogen ");
166 fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN]));
167 fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]));
168 fprintf( ioQQQ, " H+o/Hden");
169 fprintf(ioQQQ,PrintEfmt("%9.2e",hatmic ));
170 fprintf(ioQQQ,PrintEfmt("%9.2e",findspecieslocal("H-")->den/dense.gas_phase[ipHYDROGEN] ));
171 fprintf( ioQQQ, " H- H2");
172 /* this is total H2, the sum of "ground" and excited */
173 fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.H2_total/dense.gas_phase[ipHYDROGEN]));
174 fprintf(ioQQQ,PrintEfmt("%9.2e",findspecieslocal("H2+")->den/dense.gas_phase[ipHYDROGEN]));
175 fprintf( ioQQQ, " H2+ HeH+");
176 fprintf(ioQQQ,PrintEfmt("%9.2e",findspecieslocal("HeH+")->den/dense.gas_phase[ipHYDROGEN]));
177 fprintf( ioQQQ, " Ho+ ColD");
178 fprintf(ioQQQ,PrintEfmt("%9.2e",colden.colden[ipCOL_H0]));
179 fprintf(ioQQQ,PrintEfmt("%9.2e",colden.colden[ipCOL_Hp]));
180 fprintf( ioQQQ, "\n");
181
182 /* print departure coef if desired */
183 if( iso_sp[ipH_LIKE][ipHYDROGEN].lgPrtDepartCoef )
184 {
185 fprintf( ioQQQ, " Hydrogen " );
186 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].DepartCoef));
187 fprintf(ioQQQ,PrintEfmt("%9.2e", 1.));
188 fprintf( ioQQQ, " H+o/Hden");
189 fprintf(ioQQQ,PrintEfmt("%9.2e", (dense.xIonDense[ipHYDROGEN][0] + dense.xIonDense[ipHYDROGEN][1])/dense.gas_phase[ipHYDROGEN]));
190 fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.hmidep));
191 fprintf( ioQQQ, " H- H2");
192 fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.h2dep));
193 fprintf( ioQQQ, " H2+");
194 fprintf(ioQQQ,PrintEfmt("%9.2e", hmi.h2pdep));
195 fprintf( ioQQQ, " H3+");
196 fprintf(ioQQQ,PrintEfmt("%9.2e",hmi.h3pdep));
197 fprintf( ioQQQ, "\n" );
198 }
199
200 fixit(); // add a command line option to activate this
201 // print departure coefficients for diatoms species
202 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
203 {
204 if( 0 )
205 (*diatom)->H2_PrtDepartCoef();
206 }
207
208 if( prt.lgPrintHeating )
209 {
210 fprintf( ioQQQ, " ");
211 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][0]/thermal.htot));
212 fprintf( ioQQQ," ");
213 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][15]/thermal.htot));
214 fprintf( ioQQQ," ");
215 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][16]/thermal.htot));
216 fprintf( ioQQQ,"\n");
217 }
218
219 /* convert energy flux [erg cm-2 s-1] in attenuated incident continuum,
220 * diffuse emission, into equivalent temperature */
221 double TconTot = pow((rfield.EnergyIncidCont+rfield.EnergyDiffCont)/SPEEDLIGHT/7.56464e-15,0.25);
222 double Tincid = pow(rfield.EnergyIncidCont/SPEEDLIGHT/7.56464e-15,0.25);
223 double Tdiff = pow(rfield.EnergyDiffCont/SPEEDLIGHT/7.56464e-15,0.25);
224
225 /* convert sum of flux into energy density, then equivalent pressure */
226 double Pcon_nT = (Tincid + Tdiff) / SPEEDLIGHT;
227 Pcon_nT /= BOLTZMANN;
228
229 if( prt.lgPrintHeating )
230 {
231 fprintf( ioQQQ, " ");
232 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][1]/thermal.htot ));
233 fprintf( ioQQQ, " ");
234 fprintf(ioQQQ,PrintEfmt("%9.2e", 0. ));
235 fprintf( ioQQQ, " BoundCom");
236 fprintf(ioQQQ,PrintEfmt("%9.2e", ionbal.CompRecoilHeatLocal/ thermal.htot));
237 fprintf( ioQQQ, " Extra:");
238 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][20]/thermal.htot));
239 fprintf( ioQQQ, " Pairs:");
240 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][21]/ thermal.htot ));
241 fprintf( ioQQQ," H-lines\n");
242 }
243
244 /* Helium */
245 if( dense.lgElmtOn[ipHELIUM] )
246 {
247 fprintf( ioQQQ, " Helium " );
248 for( i=0; i < 3; i++ )
249 {
250 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipHELIUM][i]/dense.gas_phase[ipHELIUM]) );
251 }
252
253 fprintf( ioQQQ, " HeI 2s3S");
254 fprintf(ioQQQ,PrintEfmt("%9.2e",
255 iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2s3S].Pop()/dense.gas_phase[ipHELIUM] ));
256 fprintf( ioQQQ, " Comp H,C");
257 fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmheat ));
258 fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmcool*phycon.te));
259 fprintf( ioQQQ , " Fill Fac");
260 fprintf(ioQQQ,PrintEfmt("%9.2e", geometry.FillFac));
261 fprintf( ioQQQ , " Gam1/tot");
262 fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.H_ion_frac_photo));
263 fprintf( ioQQQ, "\n");
264
265 /* option to print departure coef */
266 if( iso_sp[ipH_LIKE][ipHELIUM].lgPrtDepartCoef )
267 {
268 fprintf( ioQQQ, " Helium " );
269 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].fb[0].DepartCoef));
270 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipH_LIKE][ipHELIUM].fb[ipH1s].DepartCoef));
271 fprintf(ioQQQ,PrintEfmt("%9.2e", 1.));
272
273 fprintf( ioQQQ, " Comp H,C");
274 fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmheat ));
275 fprintf(ioQQQ,PrintEfmt("%9.2e", rfield.cmcool*phycon.te ));
276 fprintf( ioQQQ , " Fill Fac");
277 fprintf(ioQQQ,PrintEfmt("%9.2e", geometry.FillFac ));
278 fprintf( ioQQQ , " Gam1/tot");
279 fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.H_ion_frac_photo));
280 fprintf( ioQQQ, "\n");
281 }
282
283 /* print heating from He (and others) if desired
284 * entry "lines" is induced line heating
285 * 1,12 ffheat: 2,3 he triplets, 1,20 compton */
286 if( prt.lgPrintHeating )
287 {
288 /*fprintf( ioQQQ, " %10.3e%10.3e Lines:%10.2e%10.2e Compton:%10.3e FF Heatig%10.3e\n",
289 thermal.heating[1][0]/thermal.htot, thermal.heating[1][1]/
290 thermal.htot, thermal.heating[0][22]/thermal.htot, thermal.heating[1][2]/
291 thermal.htot, thermal.heating[0][19]/thermal.htot, thermal.heating[0][11]/
292 thermal.htot );*/
293 fprintf( ioQQQ, " ");
294 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[1][0]/thermal.htot));
295 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[1][1]/thermal.htot));
296 fprintf( ioQQQ, " Lines:");
297 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][22]/thermal.htot));
298 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[1][2]/thermal.htot));
299 fprintf( ioQQQ, " Compton:");
300 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][19]/thermal.htot));
301 fprintf( ioQQQ, " FFHeatig");
302 fprintf(ioQQQ,PrintEfmt("%9.2e",thermal.heating[0][11]/thermal.htot));
303 fprintf( ioQQQ, "\n");
304 }
305
306 if( dense.lgElmtOn[ipHELIUM] )
307 {
308 /* helium singlets and triplets relative to total helium gas phase density */
309 double fac = 1./dense.gas_phase[ipHELIUM];
310 fprintf( ioQQQ, " He singlet n " );
311 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe1s1S].Pop()*fac ));
312 /* singlet n=2 complex */
313 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2s1S].Pop()*fac ));
314 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2p1P].Pop()*fac ));
315 /* singlet n=3 complex */
316 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3s1S].Pop()*fac ));
317 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3p1P].Pop()*fac ));
318 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3d1D].Pop()*fac ));
319
320 fprintf( ioQQQ, " He tripl" );
321 /* triplet n=2 complex */
322 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2s3S].Pop()*fac ));
323 fprintf(ioQQQ,PrintEfmt("%9.2e",
324 iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2p3P0].Pop()*fac+
325 iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2p3P1].Pop()*fac+
326 iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe2p3P2].Pop()*fac ));
327 /* triplet n=3 complex */
328 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3s3S].Pop()*fac ));
329 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3p3P].Pop()*fac ));
330 fprintf(ioQQQ,PrintEfmt("%9.2e", iso_sp[ipHE_LIKE][ipHELIUM].st[ipHe3d3D].Pop()*fac ));
331 fprintf( ioQQQ, "\n" );
332 }
333 }
334
335 /* loop over iso sequences to see if any populations
336 * and/or departure coefficients need to be printed */
337 for( long ipISO = ipH_LIKE; ipISO < NISO; ipISO++ )
338 {
339 for( nelem=ipISO; nelem<LIMELM; ++nelem )
340 {
341 if( dense.lgElmtOn[nelem] )
342 {
343 if( iso_sp[ipISO][nelem].lgPrtLevelPops )
344 {
345 iso_prt_pops(ipISO, nelem, false);
346 }
347 if( iso_sp[ipISO][nelem].lgPrtDepartCoef )
348 {
349 /* true says print departure coefficients
350 * instead of populations. */
351 iso_prt_pops(ipISO, nelem, true);
352 }
353 }
354 }
355 }
356
357 /* >>chng 01 dec 08, move pressure to line before grains, after radiation properties */
358 /* gas pressure, pressure due to incident radiation field, rad accel */
359 fprintf( ioQQQ, " Pressure NgasTgas");
360 fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresGasCurr/BOLTZMANN));
361 fprintf( ioQQQ, " P(total)");
362 fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresTotlCurr));
363 fprintf( ioQQQ, " P( gas )");
364 fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.PresGasCurr));
365 fprintf( ioQQQ, " P(Radtn)");
366 fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.pres_radiation_lines_curr));
367 fprintf( ioQQQ, " Rad accl");
368 fprintf(ioQQQ,PrintEfmt("%9.2e", wind.AccelTotalOutward));
369 fprintf( ioQQQ, " ForceMul");
370 fprintf(ioQQQ,PrintEfmt("%9.2e", wind.fmul));
371 fprintf( ioQQQ, "\n" );
372
373 fprintf( ioQQQ , " Texc(La) ");
374 fprintf(ioQQQ,PrintEfmt("%9.2e", hydro.TexcLya ));
375 /* attenuated incident continuum */
376 fprintf( ioQQQ , " T(I con)");
377 fprintf(ioQQQ,PrintEfmt("%9.2e", Tincid ));
378 fprintf( ioQQQ , " T(D con)");
379 fprintf(ioQQQ,PrintEfmt("%9.2e", Tdiff ));
380 fprintf( ioQQQ , " T(U tot)");
381 fprintf(ioQQQ,PrintEfmt("%9.2e", TconTot ));
382 /* print the total radiation density expressed as an equivalent gas pressure */
383 fprintf( ioQQQ , " nT (c+d)");
384 fprintf(ioQQQ,PrintEfmt("%9.2e", Pcon_nT ));
385 /* print the radiation to gas pressure */
386 fprintf( ioQQQ , " Prad/Gas");
387 fprintf(ioQQQ,PrintEfmt("%9.2e", pressure.pbeta ));
388 /* magnetic to gas pressure ratio */
389 fprintf( ioQQQ , " Pmag/Gas");
390 fprintf(ioQQQ,PrintEfmt("%9.2e", magnetic.pressure / pressure.PresGasCurr) );
391 fprintf( ioQQQ, "\n" );
392
393 if( gv.lgGrainPhysicsOn )
394 {
395 for( size_t nd=0; nd < gv.bin.size(); nd++ )
396 {
397 /* Change things so the quantum heated dust species are marked with an
398 * asterisk just after the name (K Volk)
399 * added QHMARK here and in the write statement */
400 chQHMark = (char)(( gv.bin[nd]->lgQHeat && gv.bin[nd]->lgUseQHeat ) ? '*' : ' ');
401 fprintf( ioQQQ, "%-12.12s%c DustTemp",gv.bin[nd]->chDstLab, chQHMark);
402 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->tedust));
403 fprintf( ioQQQ, " Pot Volt");
404 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->dstpot*EVRYD));
405 fprintf( ioQQQ, " Chrg (e)");
406 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->AveDustZ));
407 fprintf( ioQQQ, " drf cm/s");
408 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->DustDftVel));
409 fprintf( ioQQQ, " Heating:");
410 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->GasHeatPhotoEl));
411 fprintf( ioQQQ, " Frac tot");
412 fprintf(ioQQQ,PrintEfmt("%9.2e", gv.bin[nd]->GasHeatPhotoEl/thermal.htot));
413 fprintf( ioQQQ, "\n" );
414 }
415 }
416 /* >>chng 00 apr 20, moved save-out of quantum heating data to qheat(), by PvH */
417
418 /* heavy element molecules */
419 if( findspecieslocal("CO")->den > 0. )
420 {
421 fprintf( ioQQQ, " Molecules CH/Ctot:");
422 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("CH")->den/dense.gas_phase[ipCARBON]));
423 fprintf( ioQQQ, " CH+/Ctot");
424 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("CH+")->den/dense.gas_phase[ipCARBON]));
425 fprintf( ioQQQ, " CO/Ctot:");
426 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("CO")->den/dense.gas_phase[ipCARBON]));
427 fprintf( ioQQQ, " CO+/Ctot");
428 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("CO+")->den/dense.gas_phase[ipCARBON]));
429 fprintf( ioQQQ, " H2O/Otot");
430 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("H2O")->den/dense.gas_phase[ipOXYGEN]));
431 fprintf( ioQQQ, " OH/Ototl");
432 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("OH")->den/dense.gas_phase[ipOXYGEN]));
433 fprintf( ioQQQ, "\n");
434 }
435
436 /* information about the large H2 molecule - this just returns if not turned on */
437 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
438 (*diatom)->H2_Prt_Zone();
439
440 /* Lithium, Beryllium */
441 if( dense.lgElmtOn[ipLITHIUM] || dense.lgElmtOn[ipBERYLLIUM] ||
442 (secondaries.csupra[ipHYDROGEN][0]>0.) )
443 {
444 fprintf( ioQQQ, " Lithium " );
445 for( i=0; i < 4; i++ )
446 {
447 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipLITHIUM][i]/MAX2(1e-35,dense.gas_phase[ipLITHIUM]) ));
448 }
449 fprintf( ioQQQ, " Berylliu" );
450 for( i=0; i < 5; i++ )
451 {
452 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipBERYLLIUM][i]/MAX2(1e-35,dense.gas_phase[ipBERYLLIUM])) );
453 }
454
455 /* print secondary ionization rate for atomic hydrogen */
456 fprintf( ioQQQ, " sec ion:" );
457 fprintf(ioQQQ,PrintEfmt("%9.2e", secondaries.csupra[ipHYDROGEN][0]) );
458 fprintf( ioQQQ, "\n" );
459
460 /* option to print heating due to these stages*/
461 if( prt.lgPrintHeating )
462 {
463 fprintf( ioQQQ, " " );
464 for( i=0; i < 3; i++ )
465 {
466 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipLITHIUM][i]/ thermal.htot) );
467 }
468 fprintf( ioQQQ, " " );
469
470 for( i=0; i < 4; i++ )
471 {
472 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipBERYLLIUM][i]/thermal.htot ));
473 }
474 fprintf( ioQQQ, "\n" );
475 }
476 }
477
478 /* Boron */
479 if( dense.lgElmtOn[ipBORON] )
480 {
481 fprintf( ioQQQ, " Boron " );
482 for( i=0; i < 6; i++ )
483 {
484 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipBORON][i]/MAX2(1e-35,dense.gas_phase[ipBORON]) ));
485 }
486 fprintf( ioQQQ, "\n" );
487
488 /* option to print heating*/
489 if( prt.lgPrintHeating )
490 {
491 fprintf( ioQQQ, " " );
492 for( i=0; i < 5; i++ )
493 {
494 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipBORON][i]/thermal.htot ));
495 }
496 fprintf( ioQQQ, "\n" );
497 }
498 }
499
500 /* Carbon */
501 fprintf( ioQQQ, " Carbon " );
502 for( i=0; i < 7; i++ )
503 {
504 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[ipCARBON][i]/SDIV(dense.gas_phase[ipCARBON])) );
505 }
506 /* some molecules trail the line */
507 fprintf( ioQQQ, " H2O+/O " );
508 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("H2O+")->den/MAX2(1e-35,dense.gas_phase[ipOXYGEN]) ));
509 fprintf( ioQQQ, " OH+/Otot" );
510 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("OH+")->den/ MAX2(1e-35,dense.gas_phase[ipOXYGEN]) ));
511 /* print extra heating, normally zero */
512 fprintf( ioQQQ, " Hex(tot)" );
513 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[0][20] ));
514 fprintf( ioQQQ, "\n" );
515
516 /* option to print heating*/
517 if( prt.lgPrintHeating )
518 {
519 fprintf( ioQQQ, " " );
520 for( i=0; i < ipCARBON+1; i++ )
521 {
522 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipCARBON][i]/ thermal.htot) );
523 }
524 fprintf( ioQQQ, "\n" );
525 }
526
527 /* Nitrogen */
528 fprintf( ioQQQ, " Nitrogen " );
529 for( i=1; i <= 8; i++ )
530 {
531 fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipNITROGEN][i-1]/ SDIV(dense.gas_phase[ipNITROGEN]) ));
532 }
533 fprintf( ioQQQ, " O2/Ototl" );
534 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("O2")->den/MAX2(1e-35,dense.gas_phase[ipOXYGEN])));
535 fprintf( ioQQQ, " O2+/Otot" );
536 fprintf(ioQQQ,PrintEfmt("%9.2e", findspecieslocal("O2+")->den/ MAX2(1e-35,dense.gas_phase[ipOXYGEN]) ));
537 fprintf( ioQQQ, "\n" );
538
539 /* option to print heating*/
540 if( prt.lgPrintHeating )
541 {
542 fprintf( ioQQQ, " " );
543 for( i=0; i < ipNITROGEN+1; i++ )
544 {
545 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipNITROGEN][i]/ thermal.htot ));
546 }
547 fprintf( ioQQQ, "\n" );
548 }
549
550# if 0
551 /* Oxygen */
552 fprintf( ioQQQ, " Oxygen " );
553 for( i=1; i <= 9; i++ )
554 {
555 fprintf(ioQQQ,PrintEfmt("%9.2e",dense.xIonDense[ipOXYGEN][i-1]/ SDIV(dense.gas_phase[ipOXYGEN]) ));
556 }
557 fprintf( ioQQQ, "\n" );
558
559 /* option to print heating*/
560 if( prt.lgPrintHeating )
561 {
562 fprintf( ioQQQ, " " );
563 for( i=0; i < ipOXYGEN+1; i++ )
564 {
565 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[ipOXYGEN][i]/ thermal.htot ));
566 }
567 fprintf( ioQQQ, "\n" );
568 }
569# endif
570
571 /* now print rest of elements inside loops */
572 /* fluorine through Magnesium */
573 for( nelem=ipOXYGEN; nelem < ipALUMINIUM; ++nelem )
574 {
575 if( dense.lgElmtOn[nelem] )
576 {
577 /* print the element name and amount of shift */
578 fprintf( ioQQQ, " %10.10s ", elementnames.chElementName[nelem]);
579
580 for( i=0; i < nelem+2; i++ )
581 {
582 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[nelem][i]/dense.gas_phase[nelem] ));
583 }
584 fprintf( ioQQQ, "\n" );
585
586 /* print heating but only if needed */
587 if( prt.lgPrintHeating )
588 {
589 fprintf( ioQQQ, " " );
590 for( i=0; i < nelem+1; i++ )
591 {
592 fprintf(ioQQQ,PrintEfmt("%9.2e", thermal.heating[nelem][i]/thermal.htot ));
593 }
594 fprintf( ioQQQ, "\n" );
595 }
596 }
597 }
598
599 /* Aluminium through Zinc */
600 for( nelem=ipALUMINIUM; nelem < LIMELM; ++nelem )
601 {
602 if( dense.lgElmtOn[nelem] )
603 {
604 /* number of ionization stages to print across the page */
605 /*@-redef@*/
606 enum {LINE=13};
607 /*@+redef@*/
608 ishift = MAX2(0,dense.IonHigh[nelem]-LINE+1);
609
610 /* print the element name and amount of shift */
611 fprintf( ioQQQ, " %10.10s%2ld ", elementnames.chElementName[nelem],ishift );
612
613 for( i=0; i < LINE; i++ )
614 {
615 fprintf(ioQQQ,PrintEfmt("%9.2e", dense.xIonDense[nelem][i+ishift]/dense.gas_phase[nelem]) );
616 }
617 fprintf( ioQQQ, "\n" );
618
619 /* print heating but only if needed */
620 if( prt.lgPrintHeating )
621 {
622 fprintf( ioQQQ, " " );
623 for( i=0; i < LINE; i++ )
624 {
625 fprintf(ioQQQ,
626 PrintEfmt("%9.2e", thermal.heating[nelem][i+ishift]/thermal.htot ));
627 }
628 fprintf( ioQQQ, "\n" );
629 }
630 }
631 }
632
633 /* call FeII print routine if large atom is turned on */
634 FeIIPrint();
635 return;
636}
637
638
t_abund abund
Definition abund.cpp:5
void FeIIPrint(void)
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define PrintEfmt(F, V)
Definition cddefines.h:1472
const int ipBERYLLIUM
Definition cddefines.h:308
const int ipNITROGEN
Definition cddefines.h:311
const int ipCARBON
Definition cddefines.h:310
const int ipALUMINIUM
Definition cddefines.h:317
const int ipOXYGEN
Definition cddefines.h:312
void PrintE93(FILE *, double)
Definition service.cpp:838
const int LIMELM
Definition cddefines.h:258
const int ipLITHIUM
Definition cddefines.h:307
const int ipBORON
Definition cddefines.h:309
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
#define MAX2
Definition cddefines.h:782
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
void fixit(void)
Definition service.cpp:991
t_colden colden
Definition colden.cpp:5
#define ipCOL_Hp
Definition colden.h:26
#define ipCOL_H0
Definition colden.h:22
t_conv conv
Definition conv.cpp:5
t_dense dense
Definition dense.cpp:24
t_dynamics dynamics
Definition dynamics.cpp:44
void DynaPrtZone(void)
t_elementnames elementnames
t_geometry geometry
Definition geometry.cpp:5
GrainVar gv
Definition grainvar.cpp:5
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_hmi hmi
Definition hmi.cpp:5
t_hydro hydro
Definition hydrogenic.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipHe1s1S
Definition iso.h:41
const int ipH1s
Definition iso.h:27
const int ipHe2s3S
Definition iso.h:44
const int ipHe3d1D
Definition iso.h:56
const int ipHe2p1P
Definition iso.h:49
const int ipHe2s1S
Definition iso.h:45
const int ipHe3p3P
Definition iso.h:54
const int ipHe3s1S
Definition iso.h:53
const int ipHE_LIKE
Definition iso.h:63
const int ipHe2p3P1
Definition iso.h:47
const int ipHe3p1P
Definition iso.h:57
const int ipHe3d3D
Definition iso.h:55
const int ipHe2p3P0
Definition iso.h:46
const int ipH_LIKE
Definition iso.h:62
const int ipHe3s3S
Definition iso.h:52
void iso_prt_pops(long ipISO, long nelem, bool lgPrtDeparCoef)
const int ipHe2p3P2
Definition iso.h:48
t_magnetic magnetic
Definition magnetic.cpp:17
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molezone * findspecieslocal(const char buf[])
ChemAtomList unresolved_atom_list
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double SPEEDLIGHT
Definition physconst.h:100
UNUSED const double BOLTZMANN
Definition physconst.h:97
UNUSED const double EVRYD
Definition physconst.h:189
t_pressure pressure
Definition pressure.cpp:5
t_prt prt
Definition prt.cpp:10
void PrtLinePres(FILE *ioPRESSURE)
void PrtZone(void)
Definition prt_zone.cpp:36
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
t_secondaries secondaries
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5
Wind wind
Definition wind.cpp:5