cloudy  trunk
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 
36 void 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:");
80  fprintf(ioQQQ," Ne:");
82  fprintf(ioQQQ," R:");
84  fprintf(ioQQQ," R-R0:");
86  fprintf(ioQQQ," dR:");
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",
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 ,
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 ");
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");
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");
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",
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)");
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 */
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_elementnames::chElementName
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
colden.h
thermal.h
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:833
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:312
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
t_dense::eden
double eden
Definition: dense.h:190
count_ptr
Definition: count_ptr.h:11
t_hmi::h2pdep
double h2pdep
Definition: hmi.h:35
dense
t_dense dense
Definition: dense.cpp:24
elementnames.h
secondaries.h
t_dynamics::lgAdvection
bool lgAdvection
Definition: dynamics.h:60
Wind::fmul
realnum fmul
Definition: wind.h:65
ipBERYLLIUM
const int ipBERYLLIUM
Definition: cddefines.h:308
rfield
t_rfield rfield
Definition: rfield.cpp:8
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
geometry.h
PrintE93
void PrintE93(FILE *, double)
Definition: service.cpp:838
ipBORON
const int ipBORON
Definition: cddefines.h:309
diatom_iter
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
conv.h
rfield.h
ipHe3s3S
const int ipHe3s3S
Definition: iso.h:52
ipLITHIUM
const int ipLITHIUM
Definition: cddefines.h:307
unresolved_atom_list
ChemAtomList unresolved_atom_list
Definition: mole_species.cpp:70
ipCARBON
const int ipCARBON
Definition: cddefines.h:310
mole.h
abund.h
ipHe3s1S
const int ipHe3s1S
Definition: iso.h:53
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
DynaPrtZone
void DynaPrtZone(void)
Definition: dynamics.cpp:2009
diatoms
vector< diatomics * > diatoms
Definition: h2.cpp:8
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:46
phycon
t_phycon phycon
Definition: phycon.cpp:6
trace.h
t_mole_global::list
MoleculeList list
Definition: mole.h:317
GrainVar::lgGrainPhysicsOn
bool lgGrainPhysicsOn
Definition: grainvar.h:479
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
t_opac::TauAbsGeo
realnum ** TauAbsGeo
Definition: opacity.h:82
dynamics.h
Wind::AccelTotalOutward
realnum AccelTotalOutward
Definition: wind.h:52
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
t_hmi::H2_total
double H2_total
Definition: hmi.h:16
iso.h
ipNITROGEN
const int ipNITROGEN
Definition: cddefines.h:311
t_abund::lgAbTaON
bool lgAbTaON
Definition: abund.h:76
Wind::lgStatic
bool lgStatic(void) const
Definition: wind.h:24
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
Wind::AccelLine
realnum AccelLine
Definition: wind.h:61
opac
t_opac opac
Definition: opacity.cpp:5
wind
Wind wind
Definition: wind.cpp:5
ipCOL_Hp
#define ipCOL_Hp
Definition: colden.h:26
nzone
long int nzone
Definition: cddefines.cpp:14
ipHe3d3D
const int ipHe3d3D
Definition: iso.h:55
radius
t_radius radius
Definition: radius.cpp:5
FeIIPrint
void FeIIPrint(void)
Definition: atom_feii.cpp:1340
ipHe3d1D
const int ipHe3d1D
Definition: iso.h:56
dense.h
mole
t_mole_local mole
Definition: mole.cpp:7
trace
t_trace trace
Definition: trace.cpp:5
prt
t_prt prt
Definition: prt.cpp:10
cddefines.h
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
t_rfield::cmcool
double cmcool
Definition: rfield.h:293
t_hmi::h3pdep
double h3pdep
Definition: hmi.h:36
ipHe3p1P
const int ipHe3p1P
Definition: iso.h:57
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:47
thermal
t_thermal thermal
Definition: thermal.cpp:5
t_trace::nTrConvg
int nTrConvg
Definition: trace.h:27
t_pressure::PresGasCurr
double PresGasCurr
Definition: pressure.h:89
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:583
t_called::lgTalk
bool lgTalk
Definition: called.h:12
Wind::AccelGravity
realnum AccelGravity
Definition: wind.h:49
t_thermal::heating
double heating[LIMELM][LIMELM]
Definition: thermal.h:158
radius.h
colden
t_colden colden
Definition: colden.cpp:5
SPEEDLIGHT
const UNUSED double SPEEDLIGHT
Definition: physconst.h:100
t_ionbal::CompRecoilHeatLocal
double CompRecoilHeatLocal
Definition: ionbal.h:152
hmi.h
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
ipALUMINIUM
const int ipALUMINIUM
Definition: cddefines.h:317
t_colden::colden
realnum colden[NCOLD]
Definition: colden.h:38
t_pressure::pbeta
realnum pbeta
Definition: pressure.h:138
MAX2
#define MAX2
Definition: cddefines.h:782
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
pressure.h
LIMELM
const int LIMELM
Definition: cddefines.h:258
t_hmi::hmidep
double hmidep
Definition: hmi.h:33
t_pressure::pres_radiation_lines_curr
double pres_radiation_lines_curr
Definition: pressure.h:101
t_rfield::cmheat
double cmheat
Definition: rfield.h:292
t_conv::nPres2Ioniz
long int nPres2Ioniz
Definition: conv.h:152
iso_prt_pops
void iso_prt_pops(long ipISO, long nelem, bool lgPrtDeparCoef)
Definition: iso_solve.cpp:386
hydro
t_hydro hydro
Definition: hydrogenic.cpp:5
t_dense::IonHigh
long int IonHigh[LIMELM+1]
Definition: dense.h:120
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
t_hydro::TexcLya
realnum TexcLya
Definition: hydrogenic.h:66
ipHe2p1P
const int ipHe2p1P
Definition: iso.h:49
abund
t_abund abund
Definition: abund.cpp:5
prt.h
hydrogenic.h
t_prt::lgPrintHeating
bool lgPrintHeating
Definition: prt.h:192
t_thermal::htot
double htot
Definition: thermal.h:149
grainvar.h
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:21
wind.h
ionbal.h
magnetic.h
t_rfield::EnergyIncidCont
realnum EnergyIncidCont
Definition: rfield.h:484
ipHe3p3P
const int ipHe3p3P
Definition: iso.h:54
hmi
t_hmi hmi
Definition: hmi.cpp:5
physconst.h
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
dynamics
t_dynamics dynamics
Definition: dynamics.cpp:44
t_radius::drad
double drad
Definition: radius.h:31
t_rfield::EnergyDiffCont
realnum EnergyDiffCont
Definition: rfield.h:485
called
t_called called
Definition: called.cpp:5
conv
t_conv conv
Definition: conv.cpp:5
t_magnetic::pressure
double pressure
Definition: magnetic.h:33
gv
GrainVar gv
Definition: grainvar.cpp:5
PrtLinePres
void PrtLinePres(FILE *ioPRESSURE)
Definition: prt_linepres.cpp:16
magnetic
t_magnetic magnetic
Definition: magnetic.cpp:17
PrintEfmt
#define PrintEfmt(F, V)
Definition: cddefines.h:1472
t_pressure::PresTotlCurr
double PresTotlCurr
Definition: pressure.h:86
fixit
void fixit(void)
Definition: service.cpp:991
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
taulines.h
t_hydro::H_ion_frac_photo
realnum H_ion_frac_photo
Definition: hydrogenic.h:83
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
t_radius::depth_mid_zone
double depth_mid_zone
Definition: radius.h:41
pressure
t_pressure pressure
Definition: pressure.cpp:5
phycon.h
geometry
t_geometry geometry
Definition: geometry.cpp:5
t_hmi::h2dep
double h2dep
Definition: hmi.h:34
ipH1s
const int ipH1s
Definition: iso.h:27
Wind::windv
realnum windv
Definition: wind.h:18
ipHe2s1S
const int ipHe2s1S
Definition: iso.h:45
t_dense::lgDenFlucOn
bool lgDenFlucOn
Definition: dense.h:244
ipCOL_H0
#define ipCOL_H0
Definition: colden.h:22
opacity.h
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
called.h
PrtZone
void PrtZone(void)
Definition: prt_zone.cpp:36
h2.h
EVRYD
const UNUSED double EVRYD
Definition: physconst.h:189
t_geometry::FillFac
realnum FillFac
Definition: geometry.h:19
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:44
BOLTZMANN
const UNUSED double BOLTZMANN
Definition: physconst.h:97
t_phycon::te
double te
Definition: phycon.h:11
NISO
const int NISO
Definition: cddefines.h:261
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:48
t_radius::Radius_mid_zone
double Radius_mid_zone
Definition: radius.h:28
mole_global
t_mole_global mole_global
Definition: mole.cpp:6
t_thermal::lgUnstable
bool lgUnstable
Definition: thermal.h:53
atomfeii.h
Wind::AccelCont
realnum AccelCont
Definition: wind.h:55
t_mole_global::num_calc
int num_calc
Definition: mole.h:314
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62