cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
save_do.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*SaveDo produce save output during calculation,
4  * chTime is 'MIDL' during calculation, 'LAST' at the end */
5 /*SaveLineStuff save optical depths or source functions for all transferred lines */
6 /*Save1Line called by SaveLineStuff to produce output for one line */
7 /*SaveLineIntensity produce the 'save lines intensity' output */
8 /* save h emission, for AGN3 chapter 4, routine is below */
9 /*SaveResults1Line do single line of output for the save results and save line intensity commands */
10 /* the number of emission lines across one line of printout */
11 /*SaveSpecial generate output for the save special command */
12 /*SaveResults save results from save results command */
13 /*SaveResults1Line do single line of output for the save results and save line intensity commands */
14 /*SaveGaunts called by save gaunts command to output gaunt factors */
15 /*FindStrongestLineLabels find strongest lines contributing to point in continuum energy mesh, output in some save commands */
16 #include "cddefines.h"
17 #include "cddrive.h"
18 #include "mean.h"
19 #include "taulines.h"
20 #include "struc.h"
21 #include "iso.h"
22 #include "hyperfine.h"
23 #include "rt.h"
24 #include "magnetic.h"
25 #include "hydrogenic.h"
26 #include "secondaries.h"
27 #include "grainvar.h"
28 #include "lines.h"
29 #include "dynamics.h"
30 #include "colden.h"
31 #include "ionbal.h"
32 #include "yield.h"
33 #include "prt.h"
34 #include "iterations.h"
35 #include "heavy.h"
36 #include "conv.h"
37 #include "geometry.h"
38 #include "called.h"
39 #include "helike.h"
40 #include "opacity.h"
41 #include "phycon.h"
42 #include "timesc.h"
43 #include "radius.h"
44 #include "monitor_results.h"
45 #include "thermal.h"
46 #include "wind.h"
47 #include "hmi.h"
48 #include "pressure.h"
49 #include "elementnames.h"
50 #include "ipoint.h"
51 #include "hcmap.h"
52 #include "input.h"
53 #include "save.h"
54 #include "warnings.h"
55 #include "grid.h"
56 #include "atmdat.h"
57 #include "h2.h"
58 #include "gammas.h"
59 #include "mole.h"
60 #include "rfield.h"
61 #include "doppvel.h"
62 #include "freebound.h"
63 #include "dense.h"
64 #include "atmdat_gaunt.h"
65 #include "generic_state.h"
66 
67 enum cont_type { CT_INCI, // incident continuum, not attenuated
68  CT_OUTW_INCI, // incident continuum, attenuated through cloud
69  CT_REFL_INCI, // incident continuum, reflected towards observer
70  CT_OUTW_DIFF, // diffuse emission, attenuated through cloud
71  CT_REFL_DIFF, // diffuse emission, reflected towards observer
72  CT_OUTW_LIN, // line emission, attenuated through cloud
73  CT_REFL_LIN, // line emission, reflected towards observer
74  CT_GRN_SIL, // diffuse emission from silicate grains
75  CT_GRN_GRA, // diffuse emission from graphite grains
76  CT_GRN_TOT }; // diffuse emission from all grains
77 
78 
79 // implements the absorption option on the
80 // set save line width command
81 inline double PrettyTransmission(long j, double transmission)
82 {
83  if( save.ResolutionAbs < realnum(0.) )
84  // option to conserve energy
85  return transmission;
86  else
87  {
88  double corr = double(save.ResolutionAbs)*rfield.widflx(j)/rfield.anu(j);
89  return max(0., 1. - (1.-transmission)*corr);
90  }
91 }
92 
93 // flxCell: calculate flux of specified continuum type at frequency anu(j)
94 inline double flxCell(long j, long nEmType, cont_type ct, bool lgForceConserve = false,
95  bool lgPrtIsotropicCont = true, const realnum *trans_coef_total = NULL)
96 {
97  double val, factor = rfield.anu2(j)*EN1RYD/rfield.widflx(j);
98  // a value < 0. indicates that energy should be conserved
99  double resfactor = ( save.Resolution < realnum(0.) || lgForceConserve ) ? factor :
100  rfield.anu(j)*double(save.Resolution)*EN1RYD;
101  switch( ct )
102  {
103  case CT_INCI:
104  val = double(rfield.flux_total_incident[nEmType][j])*factor*radius.PI4_rinner_sq;
105  break;
106  case CT_OUTW_INCI:
107  ASSERT( trans_coef_total != NULL );
108  val = flux_correct_isotropic( lgPrtIsotropicCont, nEmType, j ) * factor*radius.PI4_Radius_sq;
109  if( lgForceConserve )
110  val *= double(trans_coef_total[j]);
111  else
112  val *= PrettyTransmission( j, double(trans_coef_total[j]) );
113  break;
114  case CT_REFL_INCI:
115  val = double(rfield.ConRefIncid[nEmType][j]*geometry.covgeo)*factor*radius.PI4_rinner_sq;
116  break;
117  case CT_OUTW_DIFF:
118  val = double(rfield.ConEmitOut[nEmType][j]*geometry.covgeo)*factor*radius.PI4_Radius_sq;
119  break;
120  case CT_REFL_DIFF:
121  val = double(rfield.ConEmitReflec[nEmType][j]*geometry.covgeo)*factor*radius.PI4_rinner_sq;
122  break;
123  case CT_OUTW_LIN:
124  val = double(rfield.outlin[nEmType][j]*geometry.covgeo)*resfactor*radius.PI4_Radius_sq;
125  break;
126  case CT_REFL_LIN:
127  val = double(rfield.reflin[nEmType][j]*geometry.covgeo)*resfactor*radius.PI4_rinner_sq;
128  break;
129  case CT_GRN_SIL:
130  val = double(gv.SilicateEmission[j]*geometry.covgeo)*factor*radius.PI4_Radius_sq;
131  break;
132  case CT_GRN_GRA:
133  val = double(gv.GraphiteEmission[j]*geometry.covgeo)*factor*radius.PI4_Radius_sq;
134  break;
135  case CT_GRN_TOT:
136  val = double(gv.GrainEmission[j]*geometry.covgeo)*factor*radius.PI4_Radius_sq;
137  break;
138  default:
139  TotalInsanity();
140  }
141  return val;
142 }
143 
144 /* This routine returns the spectrum needed for Keith Arnaud's XSPEC X-Ray
145  * analysis code. It should be called after cdDrive has successfully computed a
146  * model. The calling routine must ensure that the vectors have enough space to
147  * store the resulting spectrum, given the bounds and energy resolution */
148 void cdSPEC2(
149  /* option - the type of spectrum to be returned (in photons/cm^2/s/bin)
150  *
151  * 0 the total continuum, all components outward and reflected
152  *
153  * 1 the incident continuum
154  *
155  * 2 the attenuated incident continuum
156  * 3 the reflected incident continuum
157  *
158  * 4 diffuse emission, lines + continuum, outward
159  * 5 diffuse emission, lines + continuum, reflected
160  *
161  * 6 diffuse continuous emission, outward
162  * 7 diffuse continuous emission, reflected
163  *
164  * 8 total transmitted, incident + lines and continuum
165  * 9 total reflected, incident + lines and continuum
166  *
167  *10 exp(-tau) to the illuminated face */
168  int nOption,
169 
170  /* the returned spectrum, should have rfield.nflux elements */
171  realnum ReturnedSpectrum[] )
172 {
173  DEBUG_ENTRY( "cdSPEC2()" );
174 
175  const realnum *trans_coef_total = ( nOption == 0 || nOption == 2 || nOption == 8 || nOption == 10 ) ?
176  rfield.getCoarseTransCoef() : NULL;
177 
178  for( long j = 0; j < rfield.nflux; j++ )
179  {
180  double returnval;
181  if( nOption == 0 )
182  {
183  /* the attenuated incident continuum */
184  double flxatt = flxCell(j, 0, CT_OUTW_INCI, true, true, trans_coef_total);
185 
186  /* the outward emitted continuum */
187  double conem = flxCell(j, 0, CT_OUTW_DIFF, true) +
188  flxCell(j, 0, CT_OUTW_LIN, true);
189 
190  /* the reflected continuum */
191  double flxref = flxCell(j, 0, CT_REFL_INCI, true) +
192  flxCell(j, 0, CT_REFL_DIFF, true) +
193  flxCell(j, 0, CT_REFL_LIN, true);
194 
195  returnval = flxatt + conem + flxref;
196  }
197  else if( nOption == 1 )
198  {
199  /* this is the incident continuum, col 2 of save continuum command */
200  returnval = flxCell(j, 0, CT_INCI, true);
201  }
202  else if( nOption == 2 )
203  {
204  /* the attenuated transmitted continuum, no diffuse emission,
205  * col 3 of save continuum command */
206  returnval = flxCell(j, 0, CT_OUTW_INCI, true, true, trans_coef_total);
207  }
208  else if( nOption == 3 )
209  {
210  /* reflected incident continuum, col 6 of save continuum command */
211  returnval = flxCell(j, 0, CT_REFL_INCI, true);
212  }
213  else if( nOption == 4 )
214  {
215  /* all outward diffuse emission */
216  returnval = flxCell(j, 0, CT_OUTW_DIFF, true) +
217  flxCell(j, 0, CT_OUTW_LIN, true);
218  }
219  else if( nOption == 5 )
220  {
221  /* all reflected diffuse emission */
222  returnval = flxCell(j, 0, CT_REFL_DIFF, true) +
223  flxCell(j, 0, CT_REFL_LIN, true);
224  }
225  else if( nOption == 6 )
226  {
227  /* all outward line emission */
228  returnval = flxCell(j, 0, CT_OUTW_LIN, true);
229  }
230  else if( nOption == 7 )
231  {
232  /* all reflected line emission */
233  returnval = flxCell(j, 0, CT_REFL_LIN, true);
234  }
235  else if( nOption == 8 )
236  {
237  /* total transmitted continuum */
238  returnval = flxCell(j, 0, CT_OUTW_INCI, true, true, trans_coef_total) +
239  flxCell(j, 0, CT_OUTW_DIFF, true) +
240  flxCell(j, 0, CT_OUTW_LIN, true);
241  }
242  else if( nOption == 9 )
243  {
244  /* total reflected continuum */
245  returnval = flxCell(j, 0, CT_REFL_INCI, true) +
246  flxCell(j, 0, CT_REFL_DIFF, true) +
247  flxCell(j, 0, CT_REFL_LIN, true);
248  }
249  else if( nOption == 10 )
250  {
251  /* this is exp(-tau) */
252  /* This is the TOTAL attenuation in both the continuum and lines.
253  * Jon Miller discovered that the line attenuation was missing in 07.02 */
254  ASSERT( trans_coef_total != NULL );
255  returnval = double(opac.ExpmTau[j]*trans_coef_total[j]);
256  }
257  else
258  {
259  fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
260  cdEXIT(EXIT_FAILURE);
261  }
262 
263  // convert to photons/cm^2/s/bin and normalize on inner radius
264  if( nOption != 10 )
265  returnval /= rfield.anu2(j)*EN1RYD*radius.PI4_rinner_sq/rfield.widflx(j);
266 
267  ReturnedSpectrum[j] = realnum(returnval);
268  ASSERT( ReturnedSpectrum[j] >= 0.f );
269  }
270 }
271 
272 // find strongest lines contributing to point in continuum energy mesh, output in some save commands
273 STATIC void FindStrongestLineLabels( void )
274 {
275  long low_index=0;
276  long high_index=0;
277  long j_min = 0;
278  double MaxFlux = 0.;
279  long ipMaxFlux = 0;
280  long j = 0;
281 
282  ASSERT( LineSave.ipass==1 );
283 
284  while( rfield.anumax(j_min) < RYDLAM/LineSave.lines[LineSave.SortWL[0]].wavelength() )
285  j_min++;
286 
287  for( j=0; j<rfield.nflux; j++ )
288  {
289  if( j < j_min )
290  {
291  rfield.chLineLabel[j] = " ";
292  continue;
293  }
294 
295  ASSERT( LineSave.lines[LineSave.SortWL[low_index]].wavelength() != 0. );
296 
297  while( RYDLAM/LineSave.lines[LineSave.SortWL[low_index]].wavelength() < rfield.anumin(j) && low_index < LineSave.nsum-1 )
298  {
299  low_index++;
300  if( LineSave.lines[LineSave.SortWL[low_index]].wavelength() == 0. )
301  {
302  // hit the end of real wavelengths. Pad rest of labels with spaces
303  for( long j1=j; j1<rfield.nflux; j1++ )
304  rfield.chLineLabel[j1] = " ";
305  return;
306  }
307  }
308 
309  high_index = low_index;
310  ASSERT( LineSave.lines[LineSave.SortWL[high_index]].wavelength() != 0. );
311 
312  while( RYDLAM/LineSave.lines[LineSave.SortWL[high_index]].wavelength() < rfield.anumax(j) && high_index < LineSave.nsum-1 )
313  {
314  high_index++;
315  if( LineSave.lines[LineSave.SortWL[high_index]].wavelength() == 0. )
316  {
317  high_index--;
318  break;
319  }
320  }
321  // while loop found first one greater than j bin, decrement again to get back into j bin
322  high_index--;
323 
324  ASSERT( LineSave.lines[LineSave.SortWL[low_index]].wavelength() > 0. );
325  ASSERT( LineSave.lines[LineSave.SortWL[high_index]].wavelength() > 0. );
326  ASSERT( RYDLAM/LineSave.lines[LineSave.SortWL[low_index]].wavelength() >= rfield.anumin(j) );
327  ASSERT( RYDLAM/LineSave.lines[LineSave.SortWL[high_index]].wavelength() <= rfield.anumax(j) );
328 
329  MaxFlux = 0.;
330  ipMaxFlux = 0;
331 
332  for( long k = low_index; k <= high_index; k++ )
333  {
334  size_t ipLine = LineSave.SortWL[k];
335  if( LineSave.lines[ipLine].isCollisional() ||
336  LineSave.lines[ipLine].isHeat() ||
337  LineSave.lines[ipLine].isPump() ||
338  LineSave.lines[ipLine].isNInu() ||
339  LineSave.lines[ipLine].isNFnu() ||
340  LineSave.lines[ipLine].isInwardTotal() ||
341  LineSave.lines[ipLine].isInwardContinuum() ||
342  LineSave.lines[ipLine].isInward() ||
343  LineSave.lines[ipLine].isCaseA() ||
344  LineSave.lines[ipLine].isCaseB() ||
345  LineSave.lines[ipLine].isPhoPlus() ||
346  LineSave.lines[ipLine].isPcon() ||
347  LineSave.lines[ipLine].isQH() ||
348  LineSave.lines[ipLine].isUnit() )
349  continue;
350 
351  if( LineSave.lines[ipLine].SumLine(0) > MaxFlux )
352  {
353  MaxFlux = LineSave.lines[ipLine].SumLine(0);
354  ipMaxFlux = k;
355  }
356  }
357 
358  /* line label */
359  if( ipMaxFlux > 0 )
360  rfield.chLineLabel[j] = LineSave.lines[LineSave.SortWL[ipMaxFlux]].chALab();
361  }
362 
363  return;
364 }
365 
366 // index for loop over series of SAVE commands, available across file
367 STATIC long int ipPun;
368 
370 double PrtLogLin( double value )
371 {
372  if( save.lgPrtOldStyleLogs[ipPun] )
373  return log10( SDIV(value) );
374  else
375  return value;
376 }
377 
378 /*SaveLineResults do single line of output for the save results and save line intensity commands */
379 /* the number of emission lines across one line of printout */
380 namespace
381 {
382 
383  static const int LINEWIDTH = 6;
384  class SaveLineResults
385  {
386  long ipLine;
387  const LinSv *m_lines[LINEWIDTH];
388  FILE *m_ioPUN;
389  int m_typ;
390  public:
391  void save(const LinSv *line);
392  SaveLineResults(FILE *ioPUN, int typ)
393  {
394  m_ioPUN = ioPUN;
395  ipLine = 0;
396  m_typ = typ;
397  }
398  void flush()
399  {
400  if( ipLine > 0 )
401  {
402  /* this is an option to print many emission lines across an output line,
403  * the array option, or a single column of numbers, the column option
404  * that appears on the "save results" and "save intensity" commands
405  */
406  /* usual array 6 wide */
407  for( long i=0; i < ipLine; i++ )
408  {
409  fprintf( m_ioPUN, " ");
410  m_lines[i]->prt(m_ioPUN);
411  fprintf( m_ioPUN,"\t%.3e", m_lines[i]->SumLine(m_typ) );
412  /* >>chng 02 apr 24, do not print type */
413  /* single column for input into data base */
414  if( strcmp(::save.chPunRltType,"column") == 0 )
415  fprintf( m_ioPUN, "\n" );
416  }
417  if( strcmp(::save.chPunRltType,"array ") == 0 )
418  fprintf( m_ioPUN, " \n" );
419  }
420  ipLine = 0;
421  }
422  };
423 
424  int getEmType(int ipPun)
425  {
426  DEBUG_ENTRY( "getEmType()" );
427  int nEmType = (int)save.punarg[ipPun][0];
428  ASSERT( nEmType==0 || nEmType==1 );
429 
430  if (nEmType == 1 && strncmp(rfield.chCumuType,"NONE",4) == 0)
431  {
432  fprintf(ioQQQ," Must type 'set cumulative' before using 'save cumulative' output\n");
433  cdEXIT(EXIT_FAILURE);
434  }
435 
436  return nEmType;
437  }
438 }
439 
440 /*SaveGaunts called by save gaunts command to output gaunt factors */
441 STATIC void SaveGaunts(FILE* ioPUN);
442 
443 /*SaveResults save results from save results command */
444 /*SaveResults1Line do single line of output for the save results and save line intensity commands */
445 STATIC void SaveResults(FILE* ioPUN);
446 
447 STATIC void SaveLineStuff(
448  FILE * ioPUN,
449  const char *chJob ,
450  realnum xLimit);
451 
452 /* save h emission, for chapter 4, routine is below */
453 STATIC void AGN_Hemis(FILE *ioPUN );
454 
455 /*SaveLineIntensity produce the 'save lines intensity' output */
456 STATIC void SaveLineIntensity(FILE * ioPUN , long int ipPun, realnum Threshold);
457 
458 char *chDummy;
459 
460 void SaveDo(
461  /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */
462  const char *chTime)
463 {
464  long int
465  i,
466  j;
467 
468  DEBUG_ENTRY( "SaveDo()" );
469 
470  /*
471  * the "last" option on save command, to save on last iteration,
472  * is parsed at the top of the loop in only one place.
473  * no further action is needed at all for save last to work
474  * ok throughout this routine
475  */
476 
477  /*
478  * each branch can have a test whether chTime is or is not "LAST"
479  *
480  * if( lgLastOnly ) <== print after iteration is complete
481  *
482  * if "LAST" then this is last call to routine after iteration complete
483  * save only if "LAST" when results at end of iteration are needed
484  *
485  * if( ! lgLastOnly ) <== print for every zone
486  *
487  * test for .not."LAST" is for every zone result, where you do not
488  * want to save last zone twice
489  */
490 
491  /* return if no save to do */
492  if( save.nsave < 1 )
493  {
494  return;
495  }
496 
497  bool lgLastOnly = (strcmp(chTime,"LAST") == 0);
498 
499  // sort line labels if this is last call, this avoids multiple calls if several
500  // output options need sorted labels and is safer since labels will be sorted in
501  // case new code is added that reports the strong lines. The disadvantage is that
502  // we sort even if the labels are not used
503  if( lgLastOnly )
504  {
505  // sort emission line intensities so strongest lines are reported
506  FindStrongestLineLabels();
507  }
508 
509  for( ipPun=0; ipPun < save.nsave; ipPun++ )
510  {
511  /* this global variable to remember where in the save stack we are */
512  save.ipConPun = ipPun;
513 
514  /* used to identify case where no key found */
515  bool lgNoHitFirstBranch = false;
516 
517  /* iterations.lgLastIt is true if this is last iteration
518  * lgPunLstIter set true if 'last' key occurred on save command
519  * normally is false. This will skip saving if last set and
520  * this is not last iteration */
521  /* IMPORTANT: there is a second, identical if-statement halfway
522  * down this routine. Any changes here should be copied there! */
523  const bool lgActive =
524  ( iterations.lgLastIt || !save.lgPunLstIter[ipPun] ||
525  // if the sim aborted, make sure that the punch output is still done.
526  // for MIDL output this is not very useful as all previous zones from
527  // the iteration where the abort occured will be missing, but for LAST
528  // output it is important to print at least placeholders in grid runs.
529  ( lgAbort && lgLastOnly ) ||
530  ( dynamics.lgTimeDependentStatic && dynamics.lgStatic_completed ) );
531 
532  if (lgActive)
533  {
534 
535  if( strcmp(save.chSave[ipPun],"ABUN") == 0 )
536  {
537  /* save abundances vs depth */
538  if( ! lgLastOnly )
539  {
540  fprintf( save.params[ipPun].ipPnunit, "%.2f",
541  log10(MAX2(SMALLFLOAT,dense.gas_phase[ipHYDROGEN])) );
542  for( long nelem=ipHELIUM; nelem < LIMELM; nelem++ )
543  {
544  /* >>chng 05 feb 03, protect against non-positive abundances,
545  * bug caught by Marcelo Castellanos */
546  fprintf( save.params[ipPun].ipPnunit, "\t%.2f",
547  log10(MAX2(SMALLFLOAT,dense.gas_phase[nelem])) );
548  }
549  fprintf( save.params[ipPun].ipPnunit, "\n" );
550  }
551  }
552 
553  else if( strcmp(save.chSave[ipPun],"21CM") == 0 )
554  {
555  /* save information about 21 cm line */
556  if( ! lgLastOnly )
557  {
558  fprintf( save.params[ipPun].ipPnunit,
559  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
560  /* depth, cm */
561  radius.depth_mid_zone,
562  hyperfine.Tspin21cm ,
563  phycon.te ,
564  /* temperature from Lya - 21 cm optical depth ratio */
565  3.84e-7* iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon() /
566  SDIV( HFLines[0].Emis().TauCon() ),
567  /*TexcLine( &iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ),*/
568  (*HFLines[0].Lo()).Pop() ,
569  (*HFLines[0].Hi()).Pop() ,
570  OccupationNumberLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ),
571  HFLines[0].Emis().TauCon() ,
572  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon(),
573  HFLines[0].Emis().PopOpc(),
574  /* term in () is density (cm-3) of 1s, so this is n(1s) / Ts */
575  (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop())/SDIV( hyperfine.Tspin21cm),
576  /* why was above multiplied by this following term? */
577  /* *HFLines[0].EnergyErg/BOLTZMANN/4.,*/
578  HFLines[0].Emis().TauIn(),
579  TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ) ,
580  colden.H0_ov_Tspin,
581  /*>>chng 27 mar, GS, integrated 21cm spin temperature*/
582  colden.H0_21cm_lower,
583  colden.H0_21cm_upper,
584  -HFLines[0].EnergyK() / log((colden.H0_21cm_upper/3.)/colden.H0_21cm_lower)
585  );
586  }
587  }
588 
589  else if( strcmp(save.chSave[ipPun],"AGES") == 0 )
590  {
591  /* save timescales vs depth */
592  if( ! lgLastOnly )
593  {
594  int ipCO, ipOH;
595  ipCO = findspecies("CO")->index;
596  ipOH = findspecies("OH")->index;
597  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
598  /* depth, cm */
599  radius.depth_mid_zone,
600  /* cooling timescale */
601  dense.pden*BOLTZMANN*1.5*phycon.te/ thermal.htot,
602  /* H2 destruction timescale */
603  timesc.time_H2_Dest_here,
604  /* CO destruction timescale */
605  1./SDIV((ipCO != -1) ? mole.species[ipCO].snk : 0.),
606  /* OH destruction timescale */
607  1./SDIV((ipOH != -1) ? mole.species[ipOH].snk : 0.),
608  /* H recombination timescale */
609  1./(dense.eden*2.90e-10/(phycon.te70*phycon.te10/phycon.te03)) );
610  }
611  }
612 
613  else if( strcmp(save.chSave[ipPun]," AGN") == 0 )
614  {
615  if( lgLastOnly )
616  {
617  if( strcmp( save.chSaveArgs[ipPun], "HECS" ) == 0 )
618  {
619  /* this routine is in helike.c */
620  AGN_He1_CS(save.params[ipPun].ipPnunit);
621  }
622  if( strcmp( save.chSaveArgs[ipPun], "HEMI" ) == 0 )
623  {
624  /* save h emiss, for chapter 4, routine is below */
625  AGN_Hemis(save.params[ipPun].ipPnunit);
626  }
627  else
628  {
629  fprintf( ioQQQ, " SaveDo does not recognize flag %4.4s set for AGN save. This is impossible.\n",
630  save.chSave[ipPun] );
631  ShowMe();
632  cdEXIT(EXIT_FAILURE);
633  }
634  }
635  }
636 
637  else if( strcmp(save.chSave[ipPun],"MONI") == 0 )
638  {
639  if( lgLastOnly )
640  {
641  /* save the monitor output */
642  lgCheckMonitors(save.params[ipPun].ipPnunit);
643  }
644  }
645 
646  else if( strcmp(save.chSave[ipPun],"AVER") == 0 )
647  {
648  if( lgLastOnly )
649  {
650  /* save the averages output */
651  save_average( ipPun );
652  }
653  }
654 
655  else if( strncmp(save.chSave[ipPun],"CHA",3) == 0 )
656  {
657  if( lgLastOnly )
658  {
659  /* one of the charge transfer options, all in chargtran.c */
660  ChargTranPun( save.params[ipPun].ipPnunit , save.chSave[ipPun] );
661  }
662  }
663 
664  else if( strcmp( save.chSave[ipPun],"CHIA") == 0)
665  {
666  static bool lgRunOnce = true;
667  if( lgRunOnce )
668  {
669  lgRunOnce = false;
670  // save chianti collision data in physical units
671  int ipLo = 0;
672  int ipHi = 0;
673  double fupsilon = 0.;
674  double initTemp = 3.0;
675  double finalTemp = 9.1;
676  double stepTemp = 0.2;
677  for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
678  {
679  if( dBaseSpecies[ipSpecies].database == "Chianti" )
680  {
681  fprintf(save.params[ipPun].ipPnunit,"Species\tLo\tHi\tWlAng\tAul\n");
682  for( EmissionList::iterator tr=dBaseTrans[ipSpecies].Emis().begin();
683  tr != dBaseTrans[ipSpecies].Emis().end(); ++tr)
684  {
685  ipLo = tr->Tran().ipLo();
686  ipHi = tr->Tran().ipHi();
687  fprintf( save.params[ipPun].ipPnunit,"%s\t%i\t%i\t",
688  dBaseSpecies[ipSpecies].chLabel,ipLo+1,ipHi+1);
689  fprintf( save.params[ipPun].ipPnunit,"%.5e\t%.5e",tr->Tran().WLAng() , tr->Tran().Emis().Aul() );
690  fprintf( save.params[ipPun].ipPnunit,"\n");
691  }
692  // temperature scale
693  fprintf(save.params[ipPun].ipPnunit,"Species\tLo\tHi\t");
694  for(double logtemp = initTemp;logtemp < finalTemp;logtemp = logtemp + stepTemp )
695  {
696  fprintf( save.params[ipPun].ipPnunit,"\t%2.1f",logtemp);
697  }
698  fprintf( save.params[ipPun].ipPnunit,"\n");
699  // and the collision strengths
700  for( ipHi = 1; ipHi <dBaseSpecies[ipSpecies].numLevels_max; ++ipHi )
701  {
702  for( ipLo =0; ipLo < ipHi; ++ipLo )
703  {
704  fprintf( save.params[ipPun].ipPnunit,"%s\t%i\t%i\t",
705  dBaseSpecies[ipSpecies].chLabel,ipLo+1,ipHi+1);
706  for(double logtemp = initTemp;logtemp < finalTemp;logtemp = logtemp + stepTemp )
707  {
708  fupsilon = CHIANTI_Upsilon(ipSpecies, ipELECTRON, ipHi, ipLo, exp10(logtemp));
709  fprintf( save.params[ipPun].ipPnunit,"\t%.3e",fupsilon);
710  }
711  fprintf( save.params[ipPun].ipPnunit,"\n");
712  }
713  }
714  }
715  }
716  }
717  }
718 
719  else if( strcmp(save.chSave[ipPun],"COOL") == 0 ||
720  strcmp(save.chSave[ipPun],"EACH") == 0 )
721  {
722  /* save cooling, routine in file of same name */
723  if( ! lgLastOnly )
724  CoolSave(save.params[ipPun].ipPnunit, save.chSave[ipPun]);
725  }
726 
727  else if( strcmp(save.chSave[ipPun],"DOMI") == 0 )
728  {
729  /* save dominant rates */
730  if( ! lgLastOnly )
731  {
732  molecule *debug_species = findspecies( save.chSpeciesDominantRates[ipPun].c_str() );
733  if (debug_species == null_mole)
734  {
735  fprintf( ioQQQ,"Error in SAVE DOMINANT RATES, species %s not found\n",
736  save.chSpeciesDominantRates[ipPun].c_str());
737  }
738  else
739  {
740  fprintf( save.params[ipPun].ipPnunit,
741  "%e\t%e\t", radius.depth_mid_zone, mole.species[ debug_species->index ].column );
742  vector<const molecule*> debug_list;
743  debug_list.push_back(debug_species);
744  mole_dominant_rates( debug_list, save.params[ipPun].ipPnunit,
745  true, save.nLineList[ipPun], 0.0);
746  }
747  }
748  }
749 
750  else if( strcmp( save.chSave[ipPun],"CHRT") == 0 ||
751  strcmp( save.chSave[ipPun],"CHRC") == 0 )
752  {
753  bool lgCoef = false;
754  if( strcmp( save.chSave[ipPun],"CHRC") == 0 )
755  lgCoef = true;
756 
757  /* save chemistry rates command */
758  if( ! lgLastOnly )
759  {
760  bool lgHeader, lgData;
761  if( save.lgSaveHeader(ipPun) )
762  {
763  lgHeader = true;
764  lgData = false;
765  mole_save(save.params[ipPun].ipPnunit,save.optname[ipPun].c_str(),save.chSaveArgs[ipPun],lgHeader,lgData,lgCoef,radius.depth_mid_zone);
766  save.SaveHeaderDone(ipPun);
767  }
768  lgHeader = false;
769  lgData = true;
770  mole_save(save.params[ipPun].ipPnunit,save.optname[ipPun].c_str(),save.chSaveArgs[ipPun],lgHeader,lgData,lgCoef,radius.depth_mid_zone);
771  }
772  }
773 
774  else if( strncmp(save.chSave[ipPun],"DYN" , 3) == 0 )
775  {
776  /* save dynamics xxx, information about dynamical solutions */
777  if( ! lgLastOnly )
778  DynaSave(save.params[ipPun].ipPnunit ,save.chSave[ipPun][3] );
779  }
780 
781  else if( strcmp(save.chSave[ipPun],"ENTH") == 0 )
782  {
783  if( ! lgLastOnly )
784  fprintf( save.params[ipPun].ipPnunit,
785  "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
786  radius.depth_mid_zone,
787  phycon.EnthalpyDensity,
788  phycon.EnergyExcitation,
789  phycon.EnergyIonization,
790  phycon.EnergyBinding ,
791  0.5*POW2(wind.windv)*dense.xMassDensity , /* KE */
792  5./2.*pressure.PresGasCurr , /* thermal plus PdV work */
793  magnetic.EnthalpyDensity); /* magnetic terms */
794  }
795 
796  else if( strcmp(save.chSave[ipPun],"COMP") == 0 )
797  {
798  /* Compton energy exchange coefficients */
799  if( ! lgLastOnly )
800  {
801  for( long jj=0; jj<rfield.nflux; jj = jj + save.ncSaveSkip)
802  {
803  fprintf( save.params[ipPun].ipPnunit, "%10.2e%10.2e%10.2e\n",
804  rfield.anu(jj), rfield.comup[jj]/rfield.widflx(jj),
805  rfield.comdn[jj]/rfield.widflx(jj) );
806  }
807  }
808  }
809 
810  /* save continuum commands */
811  else if( strcmp(save.chSave[ipPun],"CON ") == 0 )
812  {
813  /* this is the usual "save continuum" case */
814  /* >>chng 06 apr 03, add every option to do every zone */
815  /* if lgSaveEveryZone is true then nSaveEveryZone must be positive
816  * was init to -1 */
817  bool lgPrintThis =false;
818  if( save.lgSaveEveryZone[ipPun] )
819  {
820  /* this branch, every option is on line so want to print every n zone */
821  if( ! lgLastOnly )
822  {
823  /* not last zone - print first and intermediate cases */
824  if( nzone==1 )
825  {
826  lgPrintThis = true;
827  }
828  else if( nzone%save.nSaveEveryZone[ipPun]==0 )
829  {
830  lgPrintThis = true;
831  }
832  }
833  else
834  {
835  /* this is last zone, print only if did not trip on above */
836  if( nzone!=1 && nzone%save.nSaveEveryZone[ipPun]!=0 )
837  {
838  lgPrintThis = true;
839  }
840  }
841  }
842  else
843  {
844  /* this branch, not "every", so only print the last zone */
845  if( lgLastOnly )
846  lgPrintThis = true;
847  }
848  ASSERT( !save.lgSaveEveryZone[ipPun] || save.nSaveEveryZone[ipPun]>0 );
849  if( lgPrintThis )
850  {
851  if( save.lgSaveEveryZone[ipPun] && nzone!=1)
852  fprintf( save.params[ipPun].ipPnunit, "%s\n",
853  save.chHashString );
854 
855  /* option to also print same arrays but for cumulative arrays */
856  int nEmType = getEmType(ipPun);
857 
858  const realnum *trans_coef_total=rfield.getCoarseTransCoef();
859  for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
860  {
861  /* four continua predicted here;
862  * incident, attenuated incident, emitted,
863  * then attenuated incident + emitted, last reflected
864  * reflected continuum is stored relative to inner radius
865  * others are stored for this radius */
866 
867  /* NB this code also used in save emitted,
868  * transmitted continuum commands */
869 
870  /* the incident continuum, flux_total_incident evaluated at illuminated face */
871  double flxin = flxCell(j, nEmType, CT_INCI);
872 
873  /* the reflected line emission */
874  double flxreflin = flxCell(j, nEmType, CT_REFL_LIN);
875 
876  /* the total reflected continuum, evaluated at each zone so at outer radius */
877  double flxref = flxCell(j, nEmType, CT_REFL_INCI) +
878  flxCell(j, nEmType, CT_REFL_DIFF) + flxreflin;
879 
880  /* the attenuated incident continuum, no covering factor since prediction is for view through cloud */
881  double flxatt = flxCell(j, nEmType, CT_OUTW_INCI, false,
882  save.lgPrtIsotropicCont[ipPun], trans_coef_total);
883 
884  /* outward line emission */
885  double flxlinem = flxCell(j, nEmType, CT_OUTW_LIN);
886 
887  /* the total outward emitted continuum */
888  double conem = flxCell(j, nEmType, CT_OUTW_DIFF) + flxlinem;
889 
890  /* sum of emitted and transmitted continua */
891  double flxtrn = conem + flxatt;
892 
893  /* photon energy in appropriate energy or wavelength units */
894  fprintf( save.params[ipPun].ipPnunit,"%.5e\t", AnuUnit(rfield.anu(j)) );
895  /* incident continuum */
896  fprintf( save.params[ipPun].ipPnunit,"%.3e\t", flxin );
897  /* trans cont */
898  fprintf( save.params[ipPun].ipPnunit,"%.3e\t", flxatt );
899  /* DiffOut cont */
900  fprintf( save.params[ipPun].ipPnunit,"%.3e\t", conem );
901  /* net trans cont */
902  fprintf( save.params[ipPun].ipPnunit,"%.3e\t", flxtrn );
903  /* reflected cont */
904  fprintf( save.params[ipPun].ipPnunit,"%.3e\t", flxref );
905  /* total cont */
906  fprintf( save.params[ipPun].ipPnunit,"%.3e\t", flxref + flxtrn );
907  /* reflected lines */
908  fprintf( save.params[ipPun].ipPnunit,"%.3e\t", flxreflin );
909  /* outward lines */
910  fprintf( save.params[ipPun].ipPnunit,"%.3e\t", flxlinem );
911 
912  fprintf( save.params[ipPun].ipPnunit, "%s\t%s\t",
913  /* line label */
914  rfield.chLineLabel[j].c_str() ,
915  /* cont label*/
916  rfield.chContLabel[j].c_str() );
917  /* number of lines within that cell over cell width
918  * save raw continuum has number by itself */
919  fprintf( save.params[ipPun].ipPnunit, "%.2f\n", rfield.line_count[j]/rfield.widflx(j)*rfield.anu(j) );
920  }
921  }
922  }
923 
924  else if( strcmp(save.chSave[ipPun],"CONC") == 0 )
925  {
926  /* save incident continuum */
927  /* set pointer for possible change in units of energy in continuum
928  * AnuUnit will give anu in whatever units were set with save units */
929  if( lgLastOnly )
930  {
931  /* incident continuum */
932  for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
933  {
934  double flxin = flxCell(j, 0, CT_INCI);
935  /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */
936  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.4e\t%.4e\n",
937  AnuUnit(rfield.anu(j)), flxin, rfield.OccNumbIncidCont[j]);
938  }
939  }
940  }
941 
942  else if( strcmp(save.chSave[ipPun],"CONG") == 0 )
943  {
944  /* save emitted grain continuum in optically thin limit */
945  if( lgLastOnly )
946  {
947  /* GrainMakeDiffuse broke out emission into types
948  * according to matType */
949  for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
950  {
951  double fgra = flxCell(j, 0, CT_GRN_GRA);
952  double fsil = flxCell(j, 0, CT_GRN_SIL);
953  double ftot = flxCell(j, 0, CT_GRN_TOT);
954  /* anu is .5e format to resolve energy mesh near 1 Ryd
955  * AnuUnit gives anu in whatever units were set with units option */
956  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.3e\t%.3e\t%.3e\n",
957  AnuUnit(rfield.anu(j)) , fgra, fsil, ftot );
958  }
959  }
960  }
961 
962  else if( strcmp(save.chSave[ipPun],"CONR") == 0 )
963  {
964  /* save reflected continuum */
965  /* set pointer for possible change in units of energy in continuum
966  * AnuUnit will give anu in whatever units were set with save units */
967  if( lgLastOnly )
968  {
969  if( geometry.lgSphere )
970  {
971  fprintf( save.params[ipPun].ipPnunit, " Reflected continuum not predicted when SPHERE is set.\n" );
972  fprintf( ioQQQ ,
973  "\n\n>>>>>>>>>>>>>\n Reflected continuum not predicted when SPHERE is set.\n" );
974  cdEXIT(EXIT_FAILURE);
975  }
976 
977  for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
978  {
979  // a value < 0. indicates that energy should be conserved
980  realnum resolution = ( save.Resolution < realnum(0.) ) ?
981  rfield.anu(j)/rfield.widflx(j) : save.Resolution;
982 
983  /* reflected continuum */
984  double flxref = rfield.anu2(j)*((double)rfield.ConRefIncid[0][j]+rfield.ConEmitReflec[0][j])/
985  rfield.widflx(j)*EN1RYD;
986  /* reflected lines */
987  double fref = rfield.anu(j)*resolution*rfield.reflin[0][j]*EN1RYD;
988  double av;
989  /* ratio of reflected to incident continuum, the albedo */
990  if( rfield.flux_total_incident[0][j] > 1e-25 )
991  {
992  av = rfield.ConRefIncid[0][j]/rfield.flux_total_incident[0][j];
993  }
994  else
995  {
996  av = 0.;
997  }
998  /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */
999  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4s\n",
1000  AnuUnit(rfield.anu(j)), flxref, fref, flxref + fref,
1001  av, rfield.chContLabel[j].c_str() );
1002  }
1003  }
1004  }
1005 
1006  else if( strcmp(save.chSave[ipPun],"CNVE") == 0 )
1007  {
1008  /* the save convergence error command */
1009  if( ! lgLastOnly )
1010  {
1011  fprintf( save.params[ipPun].ipPnunit,
1012  "%.5e\t%li\t%.4e\t%.4f\t%.4e\t%.4e\t%.3f\t%.4e\t%.4e\t%.4f\n",
1013  radius.depth_mid_zone,
1014  conv.nPres2Ioniz,
1015  pressure.PresTotlCurr,
1016  pressure.PresTotlError*100.,
1017  dense.EdenTrue,
1018  dense.eden,
1019  (dense.EdenTrue - dense.eden)*100./dense.EdenTrue,
1020  thermal.htot,
1021  thermal.ctot,
1022  (thermal.htot - thermal.ctot)*100./thermal.htot );
1023  }
1024  }
1025 
1026  else if( strcmp(save.chSave[ipPun],"CONB") == 0 )
1027  {
1028  /* save continuum bins binning */
1029  /* set pointer for possible change in units of energy in continuum
1030  * AnuUnit will give anu in whatever units were set with save units */
1031  if( ! lgLastOnly )
1032  {
1033  for( j=0; j<rfield.nflux_with_check; j = j + save.ncSaveSkip)
1034  {
1035  fprintf( save.params[ipPun].ipPnunit, "%14.5e\t%14.5e\t%14.5e\n",
1036  AnuUnit(rfield.anu(j)), rfield.anu(j), rfield.widflx(j) );
1037  }
1038  }
1039  }
1040 
1041  else if( strcmp(save.chSave[ipPun],"COND") == 0 )
1042  {
1043  ASSERT( fp_equal( save.punarg[ipPun][0] , (realnum)2. ) ||
1044  fp_equal( save.punarg[ipPun][0] , (realnum)1.) );
1045 
1046  /* save diffuse continuum the local line and continuous emission */
1047  if( lgLastOnly &&
1048  fp_equal(save.punarg[ipPun][0] , (realnum)1.) )
1049  {
1050  /* this option to save diffuse emission for last zone
1051  * as series of columns */
1052  for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
1053  {
1054  // a value < 0. indicates that energy should be conserved
1055  realnum resolution = ( save.Resolution < realnum(0.) ) ?
1056  rfield.anu(j)/rfield.widflx(j) : save.Resolution;
1057 
1058  double EmisLin = resolution*EN1RYD*
1059  rfield.DiffuseLineEmission[j]*rfield.anu(j);
1060  double EmisCon = rfield.ConEmitLocal[nzone][j]*
1061  rfield.anu2(j)*EN1RYD/rfield.widflx(j);
1062  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.5e\t%.5e\t%.5e\n",
1063  AnuUnit(rfield.anu(j)),
1064  EmisCon ,
1065  EmisLin ,
1066  EmisCon+EmisLin);
1067  }
1068  }
1069  else if( ! lgLastOnly &&
1070  fp_equal(save.punarg[ipPun][0] , (realnum)2.) )
1071  {
1072  /* diffuse emission per zone, with each a long row */
1073  static bool lgMustPrintHeader=true;
1074  if( lgMustPrintHeader )
1075  {
1076  lgMustPrintHeader = false;
1077  fprintf( save.params[ipPun].ipPnunit, "%.5e",
1078  AnuUnit(rfield.anu(0)) );
1079  for( j=1; j<rfield.nflux; j = j+save.ncSaveSkip)
1080  {
1081  fprintf( save.params[ipPun].ipPnunit, "\t%.5e",
1082  AnuUnit(rfield.anu(j)) );
1083  }
1084  fprintf( save.params[ipPun].ipPnunit, "\n" );
1085  }
1086  // a value < 0. indicates that energy should be conserved
1087  realnum resolution = ( save.Resolution < realnum(0.) ) ?
1088  rfield.anu(0)/rfield.widflx(0) : save.Resolution;
1089  double EmisLin = resolution*EN1RYD*
1090  rfield.DiffuseLineEmission[0]*rfield.anu(0);
1091  double EmisCon = rfield.ConEmitLocal[nzone][0]*
1092  rfield.anu2(0)*EN1RYD/rfield.widflx(0);
1093  fprintf( save.params[ipPun].ipPnunit, "%.5e",
1094  EmisCon+EmisLin);
1095  for( j=1; j<rfield.nflux; j = j+save.ncSaveSkip)
1096  {
1097  // a value < 0. indicates that energy should be conserved
1098  resolution = ( save.Resolution < realnum(0.) ) ?
1099  rfield.anu(j)/rfield.widflx(j) : save.Resolution;
1100  double EmisLin = resolution*EN1RYD*
1101  rfield.DiffuseLineEmission[j]*rfield.anu(j);
1102  double EmisCon = rfield.ConEmitLocal[nzone][j]*
1103  rfield.anu2(j)*EN1RYD/rfield.widflx(j);
1104  fprintf( save.params[ipPun].ipPnunit, "\t%.5e",
1105  EmisCon+EmisLin);
1106  }
1107  fprintf( save.params[ipPun].ipPnunit, "\n" );
1108  }
1109  }
1110 
1111  else if( strcmp(save.chSave[ipPun],"CONE") == 0 )
1112  {
1113  /* save emitted continuum */
1114  /* set pointer for possible change in units of energy in continuum
1115  * AnuUnit will give anu in whatever units were set with save units */
1116  if( lgLastOnly )
1117  {
1118  /* save emitted continuum */
1119  for( j=0; j<rfield.nflux; j+= save.ncSaveSkip)
1120  {
1121  /* this is the reflected component */
1122  double flxref = flxCell(j, 0, CT_REFL_INCI) + flxCell(j, 0, CT_REFL_DIFF) + flxCell(j, 0, CT_REFL_LIN);
1123 
1124  /* this is the total emission in the outward direction */
1125  double conem = flxCell(j, 0, CT_OUTW_DIFF) + flxCell(j, 0, CT_OUTW_LIN);
1126 
1127  /* output: photon energy, reflected, outward, total emission
1128  * >>chng 96 oct 22, format of anu to .5e to resolve energy mesh near 1 Ryd */
1129  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.3e\t%.3e\t%.3e\t%-*.*s\t%-*.*s\n",
1130  AnuUnit(rfield.anu(j)),
1131  flxref,
1132  conem,
1133  flxref + conem,
1134  NCHLAB-1, NCHLAB-1,
1135  rfield.chLineLabel[j].c_str(),
1136  NCHLAB-1, NCHLAB-1,
1137  rfield.chContLabel[j].c_str()
1138  );
1139  }
1140  }
1141  }
1142 
1143  /* save fine continuum command */
1144  else if( strcmp(save.chSave[ipPun],"CONf") == 0 )
1145  {
1146  if( lgLastOnly )
1147  {
1148  long nu_hi , nskip;
1149  if( save.punarg[ipPun][0] > 0. )
1150  /* possible lower bounds to energy range -
1151  * 0 if not set with range option*/
1152  j = ipFineCont( save.punarg[ipPun][0] );
1153  else
1154  j = 0;
1155 
1156  /* upper limit set with range option */
1157  if( save.punarg[ipPun][1]> 0. )
1158  nu_hi = ipFineCont( save.punarg[ipPun][1]);
1159  else
1160  nu_hi = rfield.nfine;
1161 
1162  /* number of cells to bring together, default is 10 */
1163  nskip = (long)save.punarg[ipPun][2];
1164  nskip = MAX2( 1, nskip );
1165 
1166  do
1167  {
1168  realnum sum1 = rfield.fine_opt_depth[j];
1169  realnum xnu = rfield.fine_anu[j];
1170  for( long jj=1; jj<nskip; ++jj )
1171  {
1172  xnu += rfield.fine_anu[j+jj];
1173  sum1 += rfield.fine_opt_depth[j+jj];
1174  }
1175  fprintf( save.params[ipPun].ipPnunit,
1176  "%.6e\t%.3e\n",
1177  AnuUnit(xnu/nskip),
1178  sexp(sum1/nskip) );
1179  j += nskip;
1180  } while( j < nu_hi );
1181  }
1182  }
1183 
1184  else if( strcmp(save.chSave[ipPun],"CONi") == 0 )
1185  {
1186  /* save continuum interactions */
1187  /* set pointer for possible change in units of energy in continuum
1188  * AnuUnit will give anu in whatever units were set with save units */
1189 
1190  /* continuum interactions */
1191  if( ! lgLastOnly )
1192  {
1193  long i1;
1194  /* this is option to set lowest energy */
1195  if( save.punarg[ipPun][0] <= 0. )
1196  {
1197  i1 = 1;
1198  }
1199  else if( save.punarg[ipPun][0] < 100. )
1200  {
1201  i1 = ipoint(save.punarg[ipPun][0]);
1202  }
1203  else
1204  {
1205  i1 = (long int)save.punarg[ipPun][0];
1206  }
1207 
1208  double fref = 0.;
1209  double fout = 0.;
1210  double fsum = 0.;
1211  double sum = 0.;
1212  double flxin = 0.;
1213 
1214  for( j=i1-1; j < rfield.nflux; j++ )
1215  {
1216  fref += rfield.flux[0][j]*opac.opacity_abs[j];
1217  fout += rfield.otslin[j]*opac.opacity_abs[j];
1218  fsum += rfield.otscon[j]*opac.opacity_abs[j];
1219  sum += rfield.ConInterOut[j]*opac.opacity_abs[j];
1220  flxin += (rfield.outlin[0][j] + rfield.outlin_noplot[j])*opac.opacity_abs[j];
1221  }
1222  fprintf( save.params[ipPun].ipPnunit, "%10.2e\t%10.2e\t%10.2e\t%10.2e\t%10.2e\n",
1223  fref, fout, fsum, sum, flxin );
1224  }
1225  }
1226 
1227  else if( strcmp(save.chSave[ipPun],"CONI") == 0 )
1228  {
1229  /* save ionizing continuum */
1230  /* set pointer for possible change in units of energy in continuum
1231  * AnuUnit will give anu in whatever units were set with save units */
1232 
1233  if( save.lgSaveEveryZone[ipPun] || (lgLastOnly) )
1234  {
1235  /* this flag will remember whether we have ever printed anything */
1236  bool lgPrt=false;
1237  if( save.lgSaveEveryZone[ipPun] )
1238  fprintf(save.params[ipPun].ipPnunit,"#save every zone %li\n", nzone);
1239 
1240  /* save ionizing continuum command
1241  * this is option to set lowest energy,
1242  * if no number was entered then this was zero */
1243  long i1;
1244  if( save.punarg[ipPun][0] <= 0. )
1245  i1 = 1;
1246  else if( save.punarg[ipPun][0] < 100. )
1247  i1 = ipoint(save.punarg[ipPun][0]);
1248  else
1249  i1 = (long int)save.punarg[ipPun][0];
1250 
1251  double sum = 0.;
1252  for( j=i1-1; j < rfield.nflux; j++ )
1253  {
1254  double flxcor = rfield.flux[0][j] +
1255  rfield.otslin[j] +
1256  rfield.otscon[j] +
1257  rfield.ConInterOut[j] +
1258  rfield.outlin[0][j] + rfield.outlin_noplot[j];
1259 
1260  sum += flxcor*opac.opacity_abs[j];
1261  }
1262 
1263  if( sum > 0. )
1264  sum = 1./sum;
1265  else
1266  sum = 1.;
1267 
1268  double fsum = 0.;
1269 
1270  for( j=i1-1; j<rfield.nflux; ++j)
1271  {
1272  double flxcor = rfield.flux[0][j] +
1273  rfield.otslin[j] +
1274  rfield.otscon[j] +
1275  rfield.ConInterOut[j]+
1276  rfield.outlin[0][j] + rfield.outlin_noplot[j];
1277 
1278  fsum += flxcor*opac.opacity_abs[j];
1279 
1280  /* punched quantities are freq, flux, flux*cross sec,
1281  * fraction of total, integral fraction of total */
1282  double RateInter = flxcor*opac.opacity_abs[j]*sum;
1283 
1284  /* punage(ipPun,2) is lowest interaction rate to consider, def=0.01 (1 percent) */
1285  /* >>chng 01 nov 22, format to c-friendly */
1286  if( (RateInter >= save.punarg[ipPun][1]) && (flxcor > SMALLFLOAT) )
1287  {
1288  lgPrt = true;
1289  /* >>chng 96 oct 22, format of anu to 11.5 to resolve energy mesh near 1 Ryd */
1290  fprintf( save.params[ipPun].ipPnunit,
1291  "%li\t%.5e\t%.2e\t%.2e\t%.2f\t%.2f\t%.2f\t%.2f\t%.2f\t%.2e\t%.2e\t%.4s\t%.4s\n",
1292  j,
1293  AnuUnit(rfield.anu(j)),
1294  flxcor,
1295  flxcor*opac.opacity_abs[j] / rfield.widflx(j),
1296  rfield.flux[0][j]/flxcor,
1297  rfield.otslin[j]/flxcor,
1298  rfield.otscon[j]/flxcor,
1299  (rfield.outlin[0][j] + rfield.outlin_noplot[j])/flxcor,
1300  rfield.ConInterOut[j]/flxcor,
1301  RateInter,
1302  fsum*sum,
1303  rfield.chLineLabel[j].c_str(),
1304  rfield.chContLabel[j].c_str() );
1305  }
1306  }
1307  if( !lgPrt )
1308  {
1309  /* entered logical block but did not print anything */
1310  fprintf(save.params[ipPun].ipPnunit,
1311  " SaveDo, the SAVE IONIZING CONTINUUM command "
1312  "did not find a strongly interacting energy, sum and fsum were %.2e %.2e\n",
1313  sum,fsum);
1314  fprintf(save.params[ipPun].ipPnunit,
1315  " SaveDo, the low-frequency energy was %.5e Ryd\n",
1316  rfield.anu(i1-1));
1317  fprintf(save.params[ipPun].ipPnunit,
1318  " You can reset the threshold for the lowest fractional "
1319  "interaction to print with the second number of the save command\n"
1320  " The fraction was %.3f and this was too large.\n",
1321  save.punarg[ipPun][1]);
1322  }
1323  }
1324  }
1325 
1326  else if( strcmp(save.chSave[ipPun],"CONS") == 0 )
1327  {
1328  if( ! lgLastOnly )
1329  {
1330  // continuum volume emissivity and opacity as a function of radius
1331  // command was "save continuum emissivity" */
1332  if( save.ipEmisFreq[ipPun] < 0 )
1333  save.ipEmisFreq[ipPun] = ipoint(save.emisfreq[ipPun].Ryd());
1334  j = save.ipEmisFreq[ipPun]-1;
1335 
1336  fprintf( save.params[ipPun].ipPnunit,
1337  "%.14e\t%.14e\t%.5e\t%.5e\t%.5e\n",
1338  radius.Radius_mid_zone,
1339  radius.depth_mid_zone,
1340  rfield.anu2(j)*rfield.ConEmitLocal[nzone][j]/rfield.widflx(j)*EN1RYD,
1341  opac.opacity_abs[j],
1342  opac.opacity_sct[j] );
1343  }
1344  }
1345 
1346  else if( strcmp(save.chSave[ipPun],"CORA") == 0 )
1347  {
1348  /* save raw continuum */
1349  /* set pointer for possible change in units of energy in continuum
1350  * AnuUnit will give anu in whatever units were set with save units */
1351 
1352  if( lgLastOnly )
1353  {
1354  /* this option to save all raw ionizing continuum */
1355  for( j=0;j<rfield.nflux;j = j + save.ncSaveSkip)
1356  {
1357  fprintf( save.params[ipPun].ipPnunit,
1358  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\t",
1359  AnuUnit(rfield.anu(j)),
1360  rfield.flux[0][j],
1361  rfield.otslin[j],
1362  rfield.otscon[j],
1363  rfield.ConRefIncid[0][j],
1364  rfield.ConEmitReflec[0][j],
1365  rfield.ConInterOut[j],
1366  rfield.outlin[0][j]+rfield.outlin_noplot[j],
1367  rfield.ConEmitOut[0][j],
1368  rfield.chLineLabel[j].c_str(),
1369  rfield.chContLabel[j].c_str()
1370  );
1371  /* number of lines within that cell */
1372  fprintf( save.params[ipPun].ipPnunit, "%li\n", rfield.line_count[j] );
1373  }
1374  }
1375  }
1376 
1377  else if( strcmp(save.chSave[ipPun],"CONT") == 0 )
1378  {
1379  /* save transmitted continuum - this is not the main "save continuum"
1380  * command - search on "CON " above
1381  * set pointer for possible change in units of energy in continuum
1382  * AnuUnit will give anu in whatever units were set with save units */
1383 
1384  if( lgLastOnly )
1385  {
1386  fprintf( save.params[ipPun].ipPnunit, "#\n" );
1387  fprintf( save.params[ipPun].ipPnunit, "%32ld # file format version number\n",
1388  VERSION_TRNCON );
1389  fprintf( save.params[ipPun].ipPnunit, "%s # check 1\n",
1390  rfield.mesh_md5sum().c_str() );
1391  union {
1392  double x;
1393  uint32 i[2];
1394  } u;
1395  u.x = rfield.emm();
1396  if( cpu.i().big_endian() )
1397  fprintf( save.params[ipPun].ipPnunit, "%23.8x %8.8x # check 2\n",
1398  u.i[0], u.i[1] );
1399  else
1400  fprintf( save.params[ipPun].ipPnunit, "%23.8x %8.8x # check 2\n",
1401  u.i[1], u.i[0] );
1402  u.x = rfield.egamry();
1403  if( cpu.i().big_endian() )
1404  fprintf( save.params[ipPun].ipPnunit, "%23.8x %8.8x # check 3\n",
1405  u.i[0], u.i[1] );
1406  else
1407  fprintf( save.params[ipPun].ipPnunit, "%23.8x %8.8x # check 3\n",
1408  u.i[1], u.i[0] );
1409  u.x = rfield.getResolutionScaleFactor();
1410  if( cpu.i().big_endian() )
1411  fprintf( save.params[ipPun].ipPnunit, "%23.8x %8.8x # check 4\n",
1412  u.i[0], u.i[1] );
1413  else
1414  fprintf( save.params[ipPun].ipPnunit, "%23.8x %8.8x # check 4\n",
1415  u.i[1], u.i[0] );
1416  fprintf( save.params[ipPun].ipPnunit, "%32.16e # radius, -1 if not set\n",
1417  radius.lgRadiusKnown ? radius.Radius : -1. );
1418  fprintf( save.params[ipPun].ipPnunit, "%32ld # nflux\n",
1419  (rfield.nflux+save.ncSaveSkip-1)/save.ncSaveSkip );
1420  fprintf( save.params[ipPun].ipPnunit, "#\n" );
1421 
1422  const realnum *trans_coef_total = rfield.getCoarseTransCoef();
1423 
1424  /* this option to save transmitted continuum */
1425  for( j=0; j < rfield.nflux; j += save.ncSaveSkip )
1426  {
1427  /* attenuated incident continuum
1428  * >>chng 97 jul 10, remove SaveLWidth from this one only since
1429  * we must conserve energy even in lines
1430  * >>chng 07 apr 26 include transmission coefficient */
1431  double flxatt = flxCell(j, 0, CT_OUTW_INCI, true, save.lgPrtIsotropicCont[ipPun],
1432  trans_coef_total);
1433 
1434  /* >>chng 00 jan 03, above did not include all contributors.
1435  * Pasted in below from usual
1436  * save continuum command */
1437  /* >>chng 04 jul 15, removed factor of save.SaveLWidth -
1438  * this should not be there to conserve energy, as explained in hazy
1439  * where command was documented, and in comment above. caught by PvH */
1440  /* >>chng 04 jul 23, incorrect use of outlin - before multiplied by an2,
1441  * quantity should be photons per Ryd, since init quantity is
1442  * photons per cell. Must div by widflx. caught by PvH */
1443  double conem = flxCell(j, 0, CT_OUTW_DIFF, true) + flxCell(j, 0, CT_OUTW_LIN, true);
1444 
1445  /* always save in intensity units, even when running in luminosity mode */
1446  double intentrn = (conem + flxatt)/radius.PI4_Radius_sq;
1447 
1448  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.3e\t%.3e\n",
1449  rfield.anu(j), intentrn, trans_coef_total[j] );
1450  }
1451  }
1452  }
1453 
1454  else if( strcmp(save.chSave[ipPun],"CON2") == 0 )
1455  {
1456  /* save total two-photon continuum */
1457  if( lgLastOnly )
1458  {
1459  /* this option to save diffuse continuum */
1460  for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
1461  {
1462  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.5e\t%.5e\n",
1463  AnuUnit(rfield.anu(j)),
1464  rfield.TotDiff2Pht[j]/rfield.widflx(j) ,
1465  rfield.TotDiff2Pht[j]*rfield.anu2(j)*EN1RYD/rfield.widflx(j));
1466  }
1467  }
1468  }
1469 
1470  else if( strcmp(save.chSave[ipPun],"DUSE") == 0 )
1471  {
1472  /* save grain extinction - includes only grain opacity, not total */
1473  if( ! lgLastOnly )
1474  {
1475  fprintf( save.params[ipPun].ipPnunit, " %.5e\t",
1476  radius.depth_mid_zone );
1477 
1478  /* visual extinction of an extended source (like a PDR)*/
1479  fprintf( save.params[ipPun].ipPnunit, "%.2e\t" , rfield.extin_mag_V_extended);
1480 
1481  /* visual extinction of point source (star)*/
1482  fprintf( save.params[ipPun].ipPnunit, "%.2e\n" , rfield.extin_mag_V_point);
1483  }
1484  }
1485 
1486  else if( strcmp(save.chSave[ipPun],"DUSO") == 0 )
1487  {
1488  /* save grain cross sections per hydrogen, cm^2/H */
1489  if( lgLastOnly )
1490  {
1491  for( j=0; j < rfield.nflux; j++ )
1492  {
1493  double scat;
1494  fprintf( save.params[ipPun].ipPnunit,
1495  "%.5e\t%.2e\t%.2e\t%.2e\t",
1496  /* photon energy or wavelength */
1497  AnuUnit(rfield.anu(j)),
1498  /* total cross section per hydrogen cm^2/H, discount forward scattering */
1499  gv.dstab[j] + gv.dstsc[j],
1500  /* absorption cross section per H */
1501  gv.dstab[j],
1502  /* scatter, with forward discounted */
1503  gv.dstsc[j] );
1504  /* add together total scattering, discounting 1-g */
1505  scat = 0.;
1506  /* sum over all grain species */
1507  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1508  {
1509  scat += gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->dstAbund;
1510  }
1511  /* finally, scattering including effects of forward scattering */
1512  fprintf( save.params[ipPun].ipPnunit,
1513  "%.2e\t", scat );
1514  fprintf( save.params[ipPun].ipPnunit,
1515  "%.2e\n", gv.dstsc[j]/SDIV(gv.dstab[j] + gv.dstsc[j]) );
1516  }
1517  }
1518  }
1519 
1520  /* save grain abundance and save grain D/G ratio commands */
1521  else if( strcmp(save.chSave[ipPun],"DUSA") == 0 ||
1522  strcmp(save.chSave[ipPun],"DUSD") == 0 )
1523  {
1524  bool lgDGRatio = ( strcmp(save.chSave[ipPun],"DUSD") == 0 );
1525 
1526  /* grain abundance */
1527  if( ! lgLastOnly )
1528  {
1529  /* print grain header first if this has not yet been done */
1530  if( save.lgSaveHeader(ipPun) )
1531  {
1532  fprintf( save.params[ipPun].ipPnunit, "#Depth" );
1533  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1534  fprintf( save.params[ipPun].ipPnunit, "\t%s", gv.bin[nd]->chDstLab );
1535  fprintf( save.params[ipPun].ipPnunit, "\ttotal\n" );
1536  save.SaveHeaderDone(ipPun);
1537  }
1538  fprintf( save.params[ipPun].ipPnunit, " %.5e",
1539  radius.depth_mid_zone );
1540  /* grain abundance per bin in g/cm^3 */
1541  double total = 0.;
1542  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1543  {
1544  double abund = gv.bin[nd]->IntVol*gv.bin[nd]->dustp[0]*
1545  gv.bin[nd]->cnv_H_pCM3;
1546  if( lgDGRatio )
1547  abund /= dense.xMassDensity;
1548  fprintf( save.params[ipPun].ipPnunit, "\t%.3e", abund );
1549  total += abund;
1550  }
1551  fprintf( save.params[ipPun].ipPnunit, "\t%.3e\n", total );
1552  }
1553  }
1554 
1555  else if( strcmp(save.chSave[ipPun],"DUSP") == 0 )
1556  {
1557  /* grain potential */
1558  if( ! lgLastOnly )
1559  {
1560  /* do labels first if this is first zone */
1561  if( save.lgSaveHeader(ipPun) )
1562  {
1563  /* first print string giving grain id */
1564  fprintf( save.params[ipPun].ipPnunit, "#Depth" );
1565  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1566  fprintf( save.params[ipPun].ipPnunit, "\t%s", gv.bin[nd]->chDstLab );
1567  fprintf( save.params[ipPun].ipPnunit, "\n" );
1568  save.SaveHeaderDone(ipPun);
1569  }
1570  fprintf( save.params[ipPun].ipPnunit, " %.5e",
1571  radius.depth_mid_zone );
1572  /* grain potential in eV */
1573  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1574  fprintf( save.params[ipPun].ipPnunit, "\t%.3e", gv.bin[nd]->dstpot*EVRYD );
1575  fprintf( save.params[ipPun].ipPnunit, "\n" );
1576  }
1577  }
1578 
1579  else if( strcmp(save.chSave[ipPun],"DUSR") == 0 )
1580  {
1581  /* grain H2 formation rates */
1582  if( ! lgLastOnly )
1583  {
1584  if( save.lgSaveHeader(ipPun) )
1585  {
1586  /* first print string giving grain id */
1587  fprintf( save.params[ipPun].ipPnunit, "#Depth" );
1588  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1589  fprintf( save.params[ipPun].ipPnunit, "\t%s", gv.bin[nd]->chDstLab );
1590  fprintf( save.params[ipPun].ipPnunit, "\n" );
1591  save.SaveHeaderDone(ipPun);
1592  }
1593  fprintf( save.params[ipPun].ipPnunit, " %.5e",
1594  radius.depth_mid_zone );
1595  /* grain formation rate for H2 */
1596  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1597  fprintf( save.params[ipPun].ipPnunit, "\t%.3e", gv.bin[nd]->rate_h2_form_grains_used );
1598  fprintf( save.params[ipPun].ipPnunit, "\n" );
1599  }
1600  }
1601 
1602  else if( strcmp(save.chSave[ipPun],"DUST") == 0 )
1603  {
1604  /* grain temperatures - K*/
1605  if( ! lgLastOnly )
1606  {
1607  /* do labels first if this is first zone */
1608  if( save.lgSaveHeader(ipPun) )
1609  {
1610  /* first print string giving grain id */
1611  fprintf( save.params[ipPun].ipPnunit, "#Depth" );
1612  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1613  fprintf( save.params[ipPun].ipPnunit, "\t%s", gv.bin[nd]->chDstLab );
1614  fprintf( save.params[ipPun].ipPnunit, "\n" );
1615  save.SaveHeaderDone(ipPun);
1616  }
1617  fprintf( save.params[ipPun].ipPnunit, " %.5e",
1618  radius.depth_mid_zone );
1619  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1620  fprintf( save.params[ipPun].ipPnunit, "\t%.3e", gv.bin[nd]->tedust );
1621  fprintf( save.params[ipPun].ipPnunit, "\n" );
1622  }
1623  }
1624 
1625  else if( strcmp(save.chSave[ipPun],"DUSC") == 0 )
1626  {
1627  /* save grain charge - eden from grains and
1628  * charge per grain in electrons / grain */
1629  if( ! lgLastOnly )
1630  {
1631  /* do labels first if this is first zone */
1632  if( save.lgSaveHeader(ipPun) )
1633  {
1634  /* first print string giving grain id */
1635  fprintf( save.params[ipPun].ipPnunit, "#Depth\tne(grn)" );
1636  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1637  fprintf( save.params[ipPun].ipPnunit, "\t%s", gv.bin[nd]->chDstLab );
1638  fprintf( save.params[ipPun].ipPnunit, "\n" );
1639  save.SaveHeaderDone(ipPun);
1640  }
1641 
1642  fprintf( save.params[ipPun].ipPnunit, " %.5e\t%.4e",
1643  radius.depth_mid_zone ,
1644  /* electron density contributed by grains, in e/cm^3,
1645  * positive number means grain supplied free electrons */
1646  gv.TotalEden );
1647 
1648  /* average charge per grain in electrons */
1649  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1650  {
1651  fprintf( save.params[ipPun].ipPnunit, "\t%.3e", gv.bin[nd]->AveDustZ );
1652  }
1653  fprintf( save.params[ipPun].ipPnunit, "\n" );
1654  }
1655  }
1656 
1657  else if( strcmp(save.chSave[ipPun],"DUSH") == 0 )
1658  {
1659  /* grain heating */
1660  if( ! lgLastOnly )
1661  {
1662  /* save grain charge, but do labels first if this is first zone */
1663  if( save.lgSaveHeader(ipPun) )
1664  {
1665  /* first print string giving grain id */
1666  fprintf( save.params[ipPun].ipPnunit, "#Depth" );
1667  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1668  fprintf( save.params[ipPun].ipPnunit, "\t%s", gv.bin[nd]->chDstLab );
1669  fprintf( save.params[ipPun].ipPnunit, "\n" );
1670  save.SaveHeaderDone(ipPun);
1671  }
1672  fprintf( save.params[ipPun].ipPnunit, " %.5e",
1673  radius.depth_mid_zone );
1674  /* grain heating */
1675  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1676  fprintf( save.params[ipPun].ipPnunit, "\t%.3e", gv.bin[nd]->GasHeatPhotoEl );
1677  fprintf( save.params[ipPun].ipPnunit, "\n" );
1678  }
1679  }
1680 
1681  else if( strcmp(save.chSave[ipPun],"DUSV") == 0 )
1682  {
1683  /* grain drift velocities */
1684  if( ! lgLastOnly )
1685  {
1686  /* save grain velocity, but do labels first if this is first zone */
1687  if( save.lgSaveHeader(ipPun) )
1688  {
1689  /* first print string giving grain id */
1690  fprintf( save.params[ipPun].ipPnunit, "#Depth" );
1691  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1692  fprintf( save.params[ipPun].ipPnunit, "\t%s", gv.bin[nd]->chDstLab );
1693  fprintf( save.params[ipPun].ipPnunit, "\n" );
1694  save.SaveHeaderDone(ipPun);
1695  }
1696  fprintf( save.params[ipPun].ipPnunit, " %.5e",
1697  radius.depth_mid_zone );
1698  /* grain drift velocity in km/s */
1699  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1700  fprintf( save.params[ipPun].ipPnunit, "\t%.3e", gv.bin[nd]->DustDftVel*1e-5 );
1701  fprintf( save.params[ipPun].ipPnunit, "\n" );
1702  }
1703  }
1704 
1705  /* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */
1706  else if( strcmp(save.chSave[ipPun],"DUSQ") == 0 )
1707  {
1708  /* save grain Qs */
1709  if( lgLastOnly )
1710  {
1711  if( save.lgSaveHeader(ipPun) )
1712  {
1713  /* first print string giving grain id */
1714  fprintf( save.params[ipPun].ipPnunit, "#grain nu/Ryd" );
1715  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1716  fprintf( save.params[ipPun].ipPnunit, "\tQ_abs%s\tQ_scat*(1-g)%s",
1717  gv.bin[nd]->chDstLab, gv.bin[nd]->chDstLab );
1718  fprintf( save.params[ipPun].ipPnunit, "\n" );
1719  save.SaveHeaderDone(ipPun);
1720  }
1721  for( j=0; j < rfield.nflux; j++ )
1722  {
1723  fprintf( save.params[ipPun].ipPnunit, " %.5e",
1724  rfield.anu(j) );
1725  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1726  {
1727  fprintf( save.params[ipPun].ipPnunit, "\t%.3e\t%.3e",
1728  gv.bin[nd]->dstab1[j]*4./gv.bin[nd]->IntArea,
1729  gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->asym[j]*4./gv.bin[nd]->IntArea );
1730  }
1731  fprintf( save.params[ipPun].ipPnunit, "\n" );
1732  }
1733  }
1734  }
1735 
1736  else if( strcmp(save.chSave[ipPun],"ELEM") == 0 )
1737  {
1738  if( ! lgLastOnly )
1739  {
1740  realnum renorm = 1.f;
1741 
1742  /* this is the index for the atomic number on the physical scale */
1743  /* >>chng 04 nov 23, use c scale throughout */
1744  long nelem = (long int)save.punarg[ipPun][0];
1745  ASSERT( nelem >= ipHYDROGEN );
1746 
1747  /* don't do this if element is not turned on */
1748  if( dense.lgElmtOn[nelem] )
1749  {
1750  /* >>chng 04 nov 23, add density option, leave as cm-3
1751  * default is still norm to total of that element */
1752  if( save.punarg[ipPun][1] == 0 )
1753  renorm = dense.gas_phase[nelem];
1754 
1755  fprintf( save.params[ipPun].ipPnunit, " %.5e", radius.depth_mid_zone );
1756 
1757  if( nelem==ipHYDROGEN )
1758  {
1759  for( j=0; j <= (nelem + 1); ++j)
1760  {
1761  fprintf( save.params[ipPun].ipPnunit, "\t%.2e",
1762  dense.xIonDense[nelem][j]/renorm );
1763  }
1764  /* H2 */
1765  fprintf( save.params[ipPun].ipPnunit, "\t%.2e",
1766  hmi.H2_total/renorm );
1767  }
1768  /* >>chng 04 nov 23 add C and O fine structure pops */
1769  else
1770  {
1771  vector<string>& chList = save.chSaveSpecies[ipPun];
1772  for ( size_t ic = 0; ic < chList.size(); ++ic )
1773  {
1774  vector<genericState> v = matchGeneric(chList[ic].c_str(),false);
1775  double dens = 0;
1776  for (size_t j=0; j<v.size(); ++j)
1777  dens += density(v[j]);
1778  fprintf( save.params[ipPun].ipPnunit, "\t%.2e",dens/renorm);
1779  }
1780  }
1781  fprintf( save.params[ipPun].ipPnunit, "\n" );
1782  }
1783  }
1784  }
1785 
1786  else if( strcmp(save.chSave[ipPun],"FRED") == 0 )
1787  {
1788  /* set with save Fred command, this punches some stuff from
1789  * Fred Hamann's dynamics project */
1790  if( ! lgLastOnly )
1791  {
1792  /* Fred's list */
1793  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.5e\t%.3e\t%.3e\t%.3e"
1794  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
1795  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
1796  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
1797  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
1798  //"\t%.3e\t%.3e\n",
1799  radius.Radius, radius.depth ,wind.windv/1e5,
1800  wind.dvdr,
1801  dense.gas_phase[ipHYDROGEN], dense.eden , phycon.te,
1802  wind.AccelLine , wind.AccelCont ,
1803  wind.fmul ,
1804  // acceleration in this zone due to electron scattering,
1805  // if incident SED was not attenuated
1806  pressure.pinzon_PresIntegElecThin/dense.xMassDensity/radius.drad_x_fillfac ,
1807  mean.xIonMean[0][ipHYDROGEN][0][0] , mean.xIonMean[0][ipHYDROGEN][1][0] ,
1808  mean.xIonMean[0][ipHELIUM][0][0] , mean.xIonMean[0][ipHELIUM][1][0] ,
1809  mean.xIonMean[0][ipHELIUM][2][0] ,
1810  mean.xIonMean[0][ipCARBON][1][0] , mean.xIonMean[0][ipCARBON][2][0] ,
1811  mean.xIonMean[0][ipCARBON][3][0] ,
1812  mean.xIonMean[0][ipOXYGEN][0][0] , mean.xIonMean[0][ipOXYGEN][1][0] ,
1813  mean.xIonMean[0][ipOXYGEN][2][0] , mean.xIonMean[0][ipOXYGEN][3][0] ,
1814  mean.xIonMean[0][ipOXYGEN][4][0] , mean.xIonMean[0][ipOXYGEN][5][0] ,
1815  mean.xIonMean[0][ipOXYGEN][6][0] , mean.xIonMean[0][ipOXYGEN][7][0] ,
1816  dense.xIonDense[ipHYDROGEN][0] , dense.xIonDense[ipHYDROGEN][1] ,
1817  dense.xIonDense[ipHELIUM][0] , dense.xIonDense[ipHELIUM][1] ,
1818  dense.xIonDense[ipHELIUM][2] ,
1819  dense.xIonDense[ipCARBON][1] , dense.xIonDense[ipCARBON][2] ,
1820  dense.xIonDense[ipCARBON][3] ,
1821  dense.xIonDense[ipOXYGEN][0] , dense.xIonDense[ipOXYGEN][1] ,
1822  dense.xIonDense[ipOXYGEN][2] , dense.xIonDense[ipOXYGEN][3] ,
1823  dense.xIonDense[ipOXYGEN][4] , dense.xIonDense[ipOXYGEN][5] ,
1824  dense.xIonDense[ipOXYGEN][6] , dense.xIonDense[ipOXYGEN][7] ,
1825  mean.xIonMean[0][ipMAGNESIUM][1][0] , dense.xIonDense[ipMAGNESIUM][1]);
1826  }
1827  }
1828 
1829  /* save spectra in fits format */
1830  else if( strcmp(save.chSave[ipPun],"FITS") == 0 )
1831  {
1832  if( lgLastOnly )
1833  saveFITSfile( save.params[ipPun].ipPnunit, NUM_OUTPUT_TYPES );
1834  }
1835  /* save gammas (but without element) */
1836  else if( strcmp(save.chSave[ipPun],"GAMt") == 0 )
1837  {
1838  if( ! lgLastOnly )
1839  {
1840  long ns;
1841  /* save photoionization rates, with the PUNCH GAMMAS command */
1842  for( long nelem=0; nelem < LIMELM; nelem++ )
1843  {
1844  for( long ion=0; ion <= nelem; ion++ )
1845  {
1846  for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
1847  {
1848  fprintf( save.params[ipPun].ipPnunit, "%3ld%3ld%3ld%10.2e%10.2e%10.2e",
1849  nelem+1, ion+1, ns+1,
1850  ionbal.PhotoRate_Shell[nelem][ion][ns][0],
1851  ionbal.PhotoRate_Shell[nelem][ion][ns][1] ,
1852  ionbal.PhotoRate_Shell[nelem][ion][ns][2] );
1853 
1854  for( j=0; j < t_yield::Inst().nelec_eject(nelem,ion,ns); j++ )
1855  {
1856  fprintf( save.params[ipPun].ipPnunit, "%5.2f",
1857  t_yield::Inst().elec_eject_frac(nelem,ion,ns,j) );
1858  }
1859  fprintf( save.params[ipPun].ipPnunit, "\n" );
1860  }
1861  }
1862  }
1863  }
1864  }
1865 
1866  /* save gammas element, ion */
1867  else if( strcmp(save.chSave[ipPun],"GAMe") == 0 )
1868  {
1869  if( ! lgLastOnly )
1870  {
1871  int ns;
1872  long nelem = (long)save.punarg[ipPun][0];
1873  long ion = (long)save.punarg[ipPun][1];
1874  /* valence shell */
1875  ns = Heavy.nsShells[nelem][ion]-1;
1876  /* show what some of the ionization sources are */
1877  GammaPrt(
1878  opac.ipElement[nelem][ion][ns][0] ,
1879  opac.ipElement[nelem][ion][ns][1] ,
1880  opac.ipElement[nelem][ion][ns][2] ,
1881  save.params[ipPun].ipPnunit,
1882  ionbal.PhotoRate_Shell[nelem][ion][ns][0] ,
1883  ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.1 );
1884  }
1885  }
1886 
1887  else if( strcmp(save.chSave[ipPun],"GAUN") == 0 )
1888  {
1889  /* save gaunt factors */
1890  if( ! lgLastOnly )
1891  SaveGaunts(save.params[ipPun].ipPnunit);
1892  }
1893 
1894  else if( strcmp(save.chSave[ipPun],"GRID") == 0 )
1895  {
1896  // generating the SAVE GRID output has been moved to cdPrepareExit()
1897  // to make sure that the output always records any type of failure
1898  }
1899  else
1900  {
1901  //no hit this branch, key should be in next
1902  lgNoHitFirstBranch = true;
1903  }
1904  }
1905 
1906  // hack needed for code to compile with Visual Studio
1907  // keep this identical to the if-statement further up!!
1908  if( lgActive )
1909  {
1910  if( strcmp(save.chSave[ipPun],"HISp") == 0 )
1911  {
1912  /* save pressure history of current zone */
1913  if( ! lgLastOnly )
1914  {
1915  /* note if pressure convergence failure occurred in history that follows */
1916  if( !conv.lgConvPres )
1917  {
1918  fprintf( save.params[ipPun].ipPnunit,
1919  "#PROBLEM Pressure not converged iter %li zone %li density-pressure follows:\n",
1920  iteration , nzone );
1921  }
1922  /* note if temperature convergence failure occurred in history that follows */
1923  if( !conv.lgConvTemp )
1924  {
1925  fprintf( save.params[ipPun].ipPnunit,
1926  "#PROBLEM Temperature not converged iter %li zone %li density-pressure follows:\n",
1927  iteration , nzone );
1928  }
1929  for( unsigned long k=0; k < conv.hist_pres_density.size(); ++k )
1930  {
1931  /* save history of density - pressure, with correct pressure */
1932  fprintf( save.params[ipPun].ipPnunit , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
1933  iteration,
1934  nzone,
1935  conv.hist_pres_density[k],
1936  conv.hist_pres_current[k],
1937  conv.hist_pres_error[k]);
1938  }
1939  }
1940  }
1941 
1942  else if( strcmp(save.chSave[ipPun],"HISt") == 0 )
1943  {
1944  /* save temperature history of current zone */
1945  if( ! lgLastOnly )
1946  {
1947  /* note if pressure convergence failure occurred in history that follows */
1948  if( !conv.lgConvPres )
1949  {
1950  fprintf( save.params[ipPun].ipPnunit,
1951  "#PROBLEM Pressure not converged iter %li zone %li temp heat cool follows:\n",
1952  iteration , nzone );
1953  }
1954  /* note if temperature convergence failure occurred in history that follows */
1955  if( !conv.lgConvTemp )
1956  {
1957  fprintf( save.params[ipPun].ipPnunit,
1958  "#PROBLEM Temperature not converged iter %li zone %li temp heat cool follows:\n",
1959  iteration , nzone );
1960  }
1961  for( unsigned long k=0; k < conv.hist_temp_temp.size(); ++k )
1962  {
1963  /* save history of density - pressure, with correct pressure */
1964  fprintf( save.params[ipPun].ipPnunit , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
1965  iteration,
1966  nzone,
1967  conv.hist_temp_temp[k],
1968  conv.hist_temp_heat[k],
1969  conv.hist_temp_cool[k]);
1970  }
1971  }
1972  }
1973 
1974  else if( strncmp(save.chSave[ipPun],"H2",2) == 0 )
1975  {
1976  /* all save info on large H2 molecule include H2 PDR pdr */
1977  save.whichDiatomToPrint[ipPun]->H2_PunchDo( save.params[ipPun].ipPnunit , save.chSave[ipPun] , chTime, ipPun );
1978  }
1979 
1980  else if( strcmp(save.chSave[ipPun],"HEAT") == 0 )
1981  {
1982  /* save heating, routine in file of same name */
1983  if( ! lgLastOnly )
1984  SaveHeat(save.params[ipPun].ipPnunit);
1985  }
1986 
1987  else if( strncmp(save.chSave[ipPun],"HE",2) == 0 )
1988  {
1989  /* various save helium commands */
1990  /* save helium line wavelengths */
1991  if( strcmp(save.chSave[ipPun] , "HELW") == 0 )
1992  {
1993  if( lgLastOnly )
1994  {
1995  /* save helium & he-like wavelengths, first header */
1996  fprintf( save.params[ipPun].ipPnunit,
1997  "Z\tElem\t2 1P->1 1S\t2 3P1->1 1S\t2 3P2->1 1S"
1998  "\t2 3S->1 1S\t2 3P2->2 3S\t2 3P1->2 3S\t2 3P0->2 3S" );
1999  fprintf( save.params[ipPun].ipPnunit, "\n" );
2000  for( long nelem=ipHELIUM; nelem<LIMELM; ++nelem )
2001  {
2002  /* print element name, nuclear charge */
2003  fprintf( save.params[ipPun].ipPnunit, "%li\t%s",
2004  nelem+1 , elementnames.chElementSym[nelem] );
2005  /*prt_wl print floating wavelength in Angstroms, in output format */
2006  fprintf( save.params[ipPun].ipPnunit, "\t" );
2007  prt_wl( save.params[ipPun].ipPnunit ,
2008  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).WLAng() );
2009  fprintf( save.params[ipPun].ipPnunit, "\t" );
2010  prt_wl( save.params[ipPun].ipPnunit ,
2011  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P1,ipHe1s1S).WLAng() );
2012  fprintf( save.params[ipPun].ipPnunit, "\t" );
2013  prt_wl( save.params[ipPun].ipPnunit ,
2014  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P2,ipHe1s1S).WLAng() );
2015  fprintf( save.params[ipPun].ipPnunit, "\t" );
2016  prt_wl( save.params[ipPun].ipPnunit ,
2017  iso_sp[ipHE_LIKE][nelem].trans(ipHe2s3S,ipHe1s1S).WLAng() );
2018  fprintf( save.params[ipPun].ipPnunit, "\t" );
2019  prt_wl( save.params[ipPun].ipPnunit ,
2020  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P2,ipHe2s3S).WLAng() );
2021  fprintf( save.params[ipPun].ipPnunit, "\t" );
2022  prt_wl( save.params[ipPun].ipPnunit ,
2023  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P1,ipHe2s3S).WLAng() );
2024  fprintf( save.params[ipPun].ipPnunit, "\t" );
2025  prt_wl( save.params[ipPun].ipPnunit ,
2026  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P0,ipHe2s3S).WLAng() );
2027  fprintf( save.params[ipPun].ipPnunit, "\n");
2028  }
2029  }
2030  }
2031  else
2032  TotalInsanity();
2033  }
2034 
2035  /* save hummer, results needed for Lya transport, to feed into David's routine */
2036  else if( strcmp(save.chSave[ipPun],"HUMM") == 0 )
2037  {
2038  double eps = iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Aul()/
2039  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Coll().ColUL( colliders );
2040  fprintf( save.params[ipPun].ipPnunit,
2041  " %.5e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
2042  radius.depth_mid_zone,
2043  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn(),
2044  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop(),
2045  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop(),
2046  phycon.te, iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().damp(), eps );
2047  }
2048 
2049  else if( strncmp( save.chSave[ipPun] , "HYD", 3 ) == 0 )
2050  {
2051  /* various save hydrogen commands */
2052  if( strcmp(save.chSave[ipPun],"HYDc") == 0 )
2053  {
2054  if( ! lgLastOnly )
2055  {
2056  /* save hydrogen physical conditions */
2057  fprintf( save.params[ipPun].ipPnunit,
2058  " %.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
2059  radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], dense.eden,
2060  dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
2061  dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
2062  hmi.H2_total/dense.gas_phase[ipHYDROGEN],
2063  findspecieslocal("H2+")->den/dense.gas_phase[ipHYDROGEN],
2064  findspecieslocal("H3+")->den/dense.gas_phase[ipHYDROGEN],
2065  findspecieslocal("H-")->den/dense.gas_phase[ipHYDROGEN] );
2066  }
2067  }
2068 
2069  else if( strcmp(save.chSave[ipPun],"HYDi") == 0 )
2070  {
2071  if( ! lgLastOnly )
2072  {
2073  /* save hydrogen ionization
2074  * this will be total decays to ground */
2075  double RateInter = 0.;
2076  double stage = iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ColIoniz*dense.eden*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
2077  double fref = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
2078  double fout = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
2079  /* 06 aug 28, from numLevels_max to _local. */
2080  for( long ion=ipH2s; ion < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local; ion++ )
2081  {
2082  /* this is total decays to ground */
2083  RateInter +=
2084  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Aul()*
2085  (iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Ploss());
2086  /* total photo from all levels */
2087  fref += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ion].gamnc*iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
2088  /* total col ion from all levels */
2089  stage += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ion].ColIoniz*dense.eden*
2090  iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
2091  fout += iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
2092  }
2093 
2094  /* make these relative to parent ion */
2095  stage /= dense.xIonDense[ipHYDROGEN][1];
2096  fref /= dense.xIonDense[ipHYDROGEN][1];
2097  fout /= dense.xIonDense[ipHYDROGEN][1];
2098 
2099  fprintf( save.params[ipPun].ipPnunit, "hion\t%4ld\t%.2e\t%.2e\t%.2e",
2100  nzone,
2101  /* photo and collision ion rates have units s-1 */
2102  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc,
2103  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ColIoniz* dense.EdenHCorr,
2104  ionbal.RateRecomTot[ipHYDROGEN][0] );
2105 
2106  fprintf( save.params[ipPun].ipPnunit, "\t%.2e",
2107  iso_sp[ipH_LIKE][ipHYDROGEN].RadRec_caseB );
2108 
2109  fprintf( save.params[ipPun].ipPnunit,
2110  "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
2111  dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0],
2112  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc/(ionbal.RateRecomTot[ipHYDROGEN][0]),
2113  iso_sp[ipH_LIKE][ipHYDROGEN].fb[1].RadRecomb[ipRecEsc],
2114  RateInter,
2115  fref/MAX2(1e-37,fout),
2116  stage/MAX2(1e-37,fout),
2117  /* simple H+ */
2118  safe_div( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*dense.xIonDense[ipHYDROGEN][0], dense.eden*dense.xIonDense[ipHYDROGEN][1] ),
2119  secondaries.csupra[ipHYDROGEN][0]);
2120 
2121  GammaPrt(iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon,rfield.nflux,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipOpac,
2122  save.params[ipPun].ipPnunit,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*
2123  0.05);
2124  }
2125  }
2126 
2127  else if( strcmp(save.chSave[ipPun],"HYDp") == 0 )
2128  {
2129  if( ! lgLastOnly )
2130  {
2131  /* save hydrogen populations
2132  * first give total atom and ion density [cm-3]*/
2133  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.2e\t%.2e",
2134  radius.depth_mid_zone,
2135  dense.xIonDense[ipHYDROGEN][0],
2136  dense.xIonDense[ipHYDROGEN][1] );
2137 
2138  /* next give state-specific densities [cm-3] */
2139  for( j=ipH1s; j < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local-1; j++ )
2140  {
2141  fprintf( save.params[ipPun].ipPnunit, "\t%.2e",
2142  iso_sp[ipH_LIKE][ipHYDROGEN].st[j].Pop() );
2143  }
2144  fprintf( save.params[ipPun].ipPnunit, "\n" );
2145  }
2146  }
2147 
2148  else if( strcmp(save.chSave[ipPun],"HYDl") == 0 )
2149  {
2150  if( lgLastOnly )
2151  {
2152  /* save hydrogen line
2153  * gives intensities and optical depths */
2154  for( long ipHi=1; ipHi<iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local -
2155  iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local; ++ipHi )
2156  {
2157  for( long ipLo=0; ipLo<ipHi; ++ipLo )
2158  {
2159  if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).ipCont() < 0 )
2160  continue;
2161  fprintf(save.params[ipPun].ipPnunit, "%li\t%li\t%li\t%li\t%.4e\t%.2e\n",
2162  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipHi].n(),
2163  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipHi].l(),
2164  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipLo].n(),
2165  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipLo].l(),
2166  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).EnergyRyd(),
2167  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).Emis().TauIn() );
2168  }
2169  }
2170  }
2171  }
2172 
2173  /* save hydrogen Lya - some details about Lya */
2174  else if( strcmp(save.chSave[ipPun],"HYDL") == 0 )
2175  {
2176  if( ! lgLastOnly )
2177  {
2178  /* the population ratio for Lya */
2179  double popul = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop()/SDIV(iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop());
2180  /* the excitation temperature of Lya */
2181  double texc = TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) );
2182  fprintf( save.params[ipPun].ipPnunit,
2183  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
2184  radius.depth_mid_zone,
2185  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn(),
2186  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot(),
2187  popul,
2188  texc,
2189  phycon.te,
2190  texc/phycon.te ,
2191  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pesc(),
2192  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest(),
2193  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().pump(),
2194  opac.opacity_abs[iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont()-1],
2195  opac.albedo[iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont()-1] );
2196  }
2197  }
2198 
2199  else if( strcmp(save.chSave[ipPun],"HYDr") == 0 )
2200  {
2201  /* save hydrogen recc - recombination cooling for AGN3 */
2202  TempChange(2500.f, false);
2203  while( phycon.te <= 20000. )
2204  {
2205  double r1;
2206  double ThinCoolingCaseB;
2207 
2208  r1 = HydroRecCool(1,0);
2209  ThinCoolingCaseB = exp10(((-25.859117 +
2210  0.16229407*phycon.telogn[0] +
2211  0.34912863*phycon.telogn[1] -
2212  0.10615964*phycon.telogn[2])/(1. +
2213  0.050866793*phycon.telogn[0] -
2214  0.014118924*phycon.telogn[1] +
2215  0.0044980897*phycon.telogn[2] +
2216  6.0969594e-5*phycon.telogn[3])))/phycon.te;
2217 
2218  fprintf( save.params[ipPun].ipPnunit, " %10.2e\t",
2219  phycon.te);
2220  fprintf( save.params[ipPun].ipPnunit, " %10.2e\t",
2221  (r1+ThinCoolingCaseB)/(BOLTZMANN*phycon.te) );
2222 
2223  fprintf( save.params[ipPun].ipPnunit, " %10.2e\t",
2224  r1/(BOLTZMANN*phycon.te));
2225 
2226  fprintf( save.params[ipPun].ipPnunit, " %10.2e\n",
2227  ThinCoolingCaseB/(BOLTZMANN*phycon.te));
2228 
2229  TempChange(phycon.te *2.f , false);
2230  }
2231  /* must exit since we have disturbed the solution */
2232  fprintf(ioQQQ , "save agn now exits since solution is disturbed.\n");
2233  cdEXIT( EXIT_SUCCESS );
2234  }
2235  else
2236  TotalInsanity();
2237  }
2238 
2239  else if( strcmp(save.chSave[ipPun],"IONI") == 0 )
2240  {
2241  if( lgLastOnly )
2242  {
2243  /* save mean ionization distribution */
2244  PrtMeanIon( 'i', false , save.params[ipPun].ipPnunit );
2245  }
2246  }
2247 
2248  /* save ionization rates */
2249  else if( strcmp(save.chSave[ipPun],"IONR") == 0 )
2250  {
2251  if( ! lgLastOnly )
2252  {
2253  /* this is element number */
2254  long nelem = (long)save.punarg[ipPun][0];
2255  fprintf( save.params[ipPun].ipPnunit,
2256  "%.5e\t%.4e\t%.4e",
2257  radius.depth_mid_zone,
2258  dense.eden ,
2259  dynamics.Rate);
2260  /* >>chng 04 oct 15, from nelem+2 to nelem+1 - array over read -
2261  * caught by PnH */
2262  for( long ion=0; ion<nelem+1; ++ion )
2263  {
2264  fprintf( save.params[ipPun].ipPnunit,
2265  "\t%.4e\t%.4e\t%.4e\t%.4e",
2266  dense.xIonDense[nelem][ion] ,
2267  ionbal.RateIonizTot(nelem,ion) ,
2268  ionbal.RateRecomTot[nelem][ion] ,
2269  dynamics.Source[nelem][ion] );
2270  }
2271  fprintf( save.params[ipPun].ipPnunit, "\n");
2272  }
2273  }
2274 
2275  else if( strcmp(save.chSave[ipPun]," IP ") == 0 )
2276  {
2277  if( lgLastOnly )
2278  {
2279  /* save valence shell ip's */
2280  for( long nelem=0; nelem < LIMELM; nelem++ )
2281  {
2282  int ion_big;
2283  double energy;
2284 
2285  /* this is the largest number of ion stages per line */
2286  const int NELEM_LINE = 10;
2287  /* this loop in case all ions do not fit across page */
2288  for( ion_big=0; ion_big<=nelem; ion_big += NELEM_LINE )
2289  {
2290  int ion_limit = MIN2(ion_big+NELEM_LINE-1,nelem);
2291 
2292  /* new line then element name */
2293  fprintf( save.params[ipPun].ipPnunit,
2294  "\n%2.2s", elementnames.chElementSym[nelem]);
2295 
2296  /* print ion stages across line */
2297  for( long ion=ion_big; ion <= ion_limit; ++ion )
2298  {
2299  fprintf( save.params[ipPun].ipPnunit, "\t%4ld", ion+1 );
2300  }
2301  fprintf( save.params[ipPun].ipPnunit, "\n" );
2302 
2303  /* this loop is over all shells */
2304  ASSERT( ion_limit < LIMELM );
2305  /* upper limit is number of shells in atom */
2306  for( long ips=0; ips < Heavy.nsShells[nelem][ion_big]; ips++ )
2307  {
2308 
2309  /* print shell label */
2310  fprintf( save.params[ipPun].ipPnunit, "%2.2s", Heavy.chShell[ips]);
2311 
2312  /* loop over possible ions */
2313  for( long ion=ion_big; ion<=ion_limit; ++ion )
2314  {
2315 
2316  /* does this subshell exist for this ion? break if it does not*/
2317  /*if( Heavy.nsShells[nelem][ion]<Heavy.nsShells[nelem][0] )*/
2318  if( ips >= Heavy.nsShells[nelem][ion] )
2319  break;
2320 
2321  /* array elements are shell, numb of electrons, element, 0 */
2322  energy = t_ADfA::Inst().ph1(ips,nelem-ion,nelem,0);
2323 
2324  /* now print threshold with correct format */
2325  if( energy < 10. )
2326  {
2327  fprintf( save.params[ipPun].ipPnunit, "\t%6.3f", energy );
2328  }
2329  else if( energy < 100. )
2330  {
2331  fprintf( save.params[ipPun].ipPnunit, "\t%6.2f", energy );
2332  }
2333  else if( energy < 1000. )
2334  {
2335  fprintf( save.params[ipPun].ipPnunit, "\t%6.1f", energy );
2336  }
2337  else
2338  {
2339  fprintf( save.params[ipPun].ipPnunit, "\t%6ld", (long)(energy) );
2340  }
2341  }
2342 
2343  /* put cs at end of long line */
2344  fprintf( save.params[ipPun].ipPnunit, "\n" );
2345  }
2346  }
2347  }
2348  }
2349  }
2350 
2351  else if( strcmp(save.chSave[ipPun],"LINC") == 0 )
2352  {
2353  /* save line cumulative */
2354  if( ! lgLastOnly )
2355  {
2356  save_line(save.params[ipPun].ipPnunit,"PUNC",
2357  save.lgEmergent[ipPun],ipPun);
2358  }
2359  }
2360  else if( strcmp(save.chSave[ipPun],"LINT") == 0 )
2361  {
2362  /* save line optical depth */
2363  if( ! lgLastOnly )
2364  {
2365  save_line(save.params[ipPun].ipPnunit,"PUNO",
2366  save.lgEmergent[ipPun],ipPun);
2367  }
2368  }
2369 
2370  else if( strcmp(save.chSave[ipPun],"LIND") == 0 )
2371  {
2372  /* save line data, then stop */
2373  SaveLineData(save.params[ipPun].ipPnunit);
2374  }
2375 
2376  else if( strcmp(save.chSave[ipPun],"LINL") == 0 )
2377  {
2378  /* save line labels, only run one time */
2379  static bool lgRunOnce=false;
2380  if( lgRunOnce )
2381  continue;
2382  lgRunOnce = true;
2383 
2384  bool lgPrintAll=false;
2385  /* LONG keyword on save line labels command sets this to 1 */
2386  if( save.punarg[ipPun][0]>0. )
2387  lgPrintAll = true;
2388  prt_LineLabels(save.params[ipPun].ipPnunit , lgPrintAll );
2389  }
2390 
2391  else if( strcmp(save.chSave[ipPun],"LINO") == 0 )
2392  {
2393  if( lgLastOnly )
2394  {
2395  /* save line optical depths, routine is below, file static */
2396  SaveLineStuff(save.params[ipPun].ipPnunit,"optical" , save.punarg[ipPun][0]);
2397  }
2398  }
2399 
2400  else if( strcmp(save.chSave[ipPun],"LINP") == 0 )
2401  {
2402  if( ! lgLastOnly )
2403  {
2404  static bool lgFirst=true;
2405  /* save line populations, need to do this twice if very first
2406  * call since first call to SaveLineStuff generates atomic parameters
2407  * rather than level pops, routine is below, file static */
2408  SaveLineStuff(save.params[ipPun].ipPnunit,"populat" , save.punarg[ipPun][0]);
2409  if( lgFirst )
2410  {
2411  lgFirst = false;
2412  SaveLineStuff(save.params[ipPun].ipPnunit,"populat" , save.punarg[ipPun][0]);
2413  }
2414  }
2415  }
2416 
2417  else if( strcmp(save.chSave[ipPun],"LINS") == 0 )
2418  {
2419  /* save line emissivity */
2420  if( ! lgLastOnly )
2421  {
2422  save_line(save.params[ipPun].ipPnunit,"PUNS",
2423  save.lgEmergent[ipPun],ipPun);
2424  }
2425  }
2426 
2427  else if( strcmp(save.chSave[ipPun],"LINR") == 0 )
2428  {
2429  /* save line RT */
2430  if( ! lgLastOnly )
2431  Save_Line_RT( save.params[ipPun].ipPnunit);
2432  }
2433 
2434  else if( strcmp(save.chSave[ipPun],"LINA") == 0 )
2435  {
2436  /* save line array */
2437  if( lgLastOnly )
2438  {
2439  /* save out all lines with energies */
2440  for( j=0; j < LineSave.nsum; j++ )
2441  {
2442  if( LineSave.lines[j].wavelength() > 0. &&
2443  LineSave.lines[j].SumLine(0) > 0. )
2444  {
2445  /* line energy, in units set with units option */
2446  fprintf( save.params[ipPun].ipPnunit, "%12.5e",
2447  AnuUnit((realnum)RYDLAM/LineSave.lines[j].wavelength()) );
2448  /* line label */
2449  fprintf( save.params[ipPun].ipPnunit, "\t");
2450  LineSave.lines[j].prt(save.params[ipPun].ipPnunit);
2451  /* intrinsic intensity */
2452  fprintf( save.params[ipPun].ipPnunit, "\t%8.3f",
2453  log10(SDIV(LineSave.lines[j].SumLine(0) * radius.Conv2PrtInten)) );
2454  /* emergent line intensity, r recombination */
2455  fprintf( save.params[ipPun].ipPnunit, "\t%8.3f",
2456  log10(SDIV(LineSave.lines[j].SumLine(1) * radius.Conv2PrtInten) ) );
2457  /* type of line, i for info, etc */
2458  fprintf( save.params[ipPun].ipPnunit, " \t%c\n",
2459  LineSave.lines[j].chSumTyp());
2460  }
2461  }
2462  }
2463  }
2464 
2465  else if( strcmp(save.chSave[ipPun],"LINI") == 0 )
2466  {
2467  if( lgLastOnly &&
2468  (nzone/save.LinEvery)*save.LinEvery != nzone )
2469  {
2470  /* this is the last zone
2471  * save line intensities - but do not do last zone twice */
2472  SaveLineIntensity(save.params[ipPun].ipPnunit , ipPun , save.punarg[ipPun][0] );
2473  }
2474  else if( ! lgLastOnly )
2475  {
2476  /* following so we only save first zone if LinEvery reset */
2477  if( (save.lgLinEvery && nzone == 1) ||
2478  (nzone/save.LinEvery)*save.LinEvery == nzone )
2479  {
2480  /* this is middle of calculation
2481  * save line intensities */
2482  SaveLineIntensity(save.params[ipPun].ipPnunit , ipPun , save.punarg[ipPun][0]);
2483  }
2484  }
2485  }
2486 
2487  else if( strcmp( save.chSave[ipPun],"LEIL") == 0)
2488  {
2489  /* some line intensities for the Leiden PDR,
2490  * but only do this when calculation is complete */
2491  if( lgLastOnly )
2492  {
2493  double absval , rel;
2494  long int n;
2495  /* the lines we will find,
2496  * for a sample list of PDR lines look at LineList_PDR_H2.dat
2497  * in the cloudy data dir */
2498  /* the number of H2 lines */
2499  const int NLINE_H2 = 31;
2500  /* the number of lines which are not H2 */
2501  const int NLINE_NOTH_H2 = 5;
2502  /* the labels and wavelengths for the lines that are not H2 */
2503  char chLabel[NLINE_NOTH_H2][NCHLAB]=
2504  { "C 2", "O 1", "O 1", "C 1", "C 1" };
2505  double Wl[NLINE_NOTH_H2]=
2506  { 157.636 , 63.1679 , 145.495, 609.590 , 370.269 };
2507  /* these are wavelengths in microns, conv to Angstroms before call */
2508  /* >>chng 05 sep 06, many of following wavelengths updated to agree
2509  * with output - apparently not updated when energies changed */
2510  double Wl_H2[NLINE_H2]=
2511  {2.12125,
2512  28.213 , 17.03 , 12.2752, 9.66228, 8.02362, 6.90725, 6.10718, 5.50996, 5.05148, 4.69342,
2513  4.40836, 4.17983, 3.99573, 3.84534, 3.72257, 3.62531, 3.54606, 3.48530, 3.43697, 3.40299,
2514  3.37995, 3.36794, 3.36534, 3.37087, 3.38671, 3.40989, 3.44080, 3.48530, 3.54226, 3.60346};
2515  /* print a header for the lines */
2516  for( n=0; n<NLINE_NOTH_H2; ++n )
2517  {
2518  prt_line_inlist( save.params[ipPun].ipPnunit, chLabel[n], Wl[n] );
2519  /* get the line, non positive return says didn't find it */
2520  /* arguments are 4-char label, wavelength, return log total intensity, linear rel inten */
2521  if( cdLine( chLabel[n] , (realnum)(Wl[n]*1e4) , &rel, &absval ) <= 0 )
2522  {
2523  fprintf(save.params[ipPun].ipPnunit, " did not find\n");
2524  }
2525  else
2526  {
2527  fprintf(save.params[ipPun].ipPnunit, "\t%.3e\t%.3e\n", absval, rel);
2528  }
2529  }
2530  fprintf(save.params[ipPun].ipPnunit, "\n\n\n");
2531 
2532  /* only print the H2 lines if the big molecule is turned on */
2533  if( h2.lgEnabled )
2534  {
2535  fprintf(save.params[ipPun].ipPnunit,
2536  "Here are some of the H2 Intensities, The first one is the\n"
2537  "1-0 S(0) line and the following ones are the 0-0 S(X)\n"
2538  "lines where X goes from 0 to 29\n\n");
2539  for( n=0; n<NLINE_H2; ++n )
2540  {
2541  prt_line_inlist( save.params[ipPun].ipPnunit, "H2 ", Wl_H2[n] );
2542  /* get the line, non positive return says didn't find it */
2543  if( cdLine( "H2 " , (realnum)(Wl_H2[n]*1e4) , &rel, &absval ) <= 0 )
2544  {
2545  fprintf(save.params[ipPun].ipPnunit, " did not find\n");
2546  }
2547  else
2548  {
2549  fprintf(save.params[ipPun].ipPnunit, "\t%.3e\t%.3e\n", absval, rel);
2550  }
2551  }
2552  }
2553  }
2554  }
2555 
2556  else if( strcmp( save.chSave[ipPun],"LEIS") == 0)
2557  {
2558  if( ! lgLastOnly )
2559  {
2560  /* get some column densities we shall need */
2561  double col_ci , col_oi , col_cii, col_heii;
2562  if( cdColm("carb" , 1 , &col_ci ) )
2563  TotalInsanity();
2564  if( cdColm("carb" , 2 , &col_cii ) )
2565  TotalInsanity();
2566  if( cdColm("oxyg" , 1 , &col_oi ) )
2567  TotalInsanity();
2568  if( cdColm("heli" , 2 , &col_heii ) )
2569  TotalInsanity();
2570  /* save Leiden structure - some numbers for the Leiden PDR model comparisons */
2571  fprintf( save.params[ipPun].ipPnunit,
2572  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
2573  "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
2574  "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
2575  "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
2576  "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
2577  /* depth of this point */
2578  radius.depth_mid_zone,
2579  /* A_V for an extended source */
2580  0.00,
2581  /* A_V for a point source */
2582  rfield.extin_mag_V_point,
2583  /* temperature */
2584  phycon.te ,
2585  dense.xIonDense[ipHYDROGEN][0],
2586  hmi.H2_total,
2587  dense.xIonDense[ipCARBON][0],
2588  dense.xIonDense[ipCARBON][1],
2589  dense.xIonDense[ipOXYGEN][0],
2590  findspecieslocal("CO")->den,
2591  findspecieslocal("O2")->den,
2592  findspecieslocal("CH")->den,
2593  findspecieslocal("OH")->den,
2594  dense.eden,
2595  dense.xIonDense[ipHELIUM][1],
2596  dense.xIonDense[ipHYDROGEN][1],
2597  findspecieslocal("H3+")->den,
2598  findspecieslocal("H")->column,
2599  findspecieslocal("H2")->column+
2600  findspecieslocal("H2*")->column,
2601  col_ci,
2602  col_cii,
2603  col_oi,
2604  findspecieslocal("CO")->column,
2605  findspecieslocal("O2")->column,
2606  findspecieslocal("CH")->column,
2607  findspecieslocal("OH")->column,
2608  colden.colden[ipCOL_elec],
2609  col_heii,
2610  findspecieslocal("H+")->column,
2611  findspecieslocal("H3+")->column,
2612  hmi.H2_Solomon_dissoc_rate_used_H2g ,
2613  gv.rate_h2_form_grains_used_total,
2614  hmi.H2_photodissoc_used_H2g,
2615  hmi.UV_Cont_rel2_Draine_DB96_depth,
2616  /* CO and C dissociation rate */
2617  mole.findrk("PHOTON,CO=>C,O"),
2618  /* total CI ionization rate */
2619  ionbal.PhotoRate_Shell[ipCARBON][0][2][0],
2620  /* total heating, erg cm-3 s-1 */
2621  thermal.htot,
2622  /* total cooling, erg cm-3 s-1 */
2623  thermal.ctot,
2624  /* GrnP grain photo heating */
2625  thermal.heating(0,13),
2626  /* grain collisional cooling */
2627  MAX2(0.,gv.GasCoolColl),
2628  /* grain collisional heating */
2629  -1.*MIN2(0.,gv.GasCoolColl),
2630  /* COds - CO dissociation heating */
2631  thermal.heating(0,9),
2632  /* H2dH-Heating due to H2 dissociation */
2633  hmi.HeatH2Dish_used,
2634  /* H2vH-Heating due to collisions with H2 */
2635  hmi.HeatH2Dexc_used ,
2636  /* ChaT - charge transfer heating */
2637  thermal.heating(0,24) ,
2638  /* cosmic ray heating */
2639  thermal.heating(1,6) ,
2640  /* heating due to atoms of various heavy elements */
2641  thermal.heating(ipMAGNESIUM,0),
2642  thermal.heating(ipSULPHUR,0),
2643  thermal.heating(ipSILICON,0),
2644  thermal.heating(ipIRON,0),
2645  thermal.heating(ipSODIUM,0),
2646  thermal.heating(ipALUMINIUM,0),
2647  thermal.heating(ipCARBON,0),
2648  0.0,
2649  0.0,
2650  0.0,
2651  0.0,
2652  0.0);
2653  }
2654  }
2655 
2656  else if( strcmp( save.chSave[ipPun],"LLST") == 0)
2657  {
2658  /* save linelist command - do on last iteration */
2659  if( lgLastOnly )
2660  {
2661  fprintf( save.params[ipPun].ipPnunit, "iteration %li" , iteration );
2662  if( save.punarg[ipPun][1] )// column print
2663  fprintf( save.params[ipPun].ipPnunit, "\n" );
2664 
2665  /* -1 is flag saying that this save command was not set */
2666  if( save.nLineList[ipPun] < 0 )
2667  TotalInsanity();
2668 
2669  int LineType = 0;
2670  if( save.lgEmergent[ipPun] )
2671  LineType = 1;
2672  if( save.lgCumulative[ipPun] )
2673  LineType += 2;
2674 
2675  bool lgBadLine = false;
2676  /* loop over all lines in the file we read */
2677  for( j=0; j<save.nLineList[ipPun]; ++j )
2678  {
2679  double relative , absolute, PrtQuantity;
2680  if( (cdLine( save.chLineListLabel[ipPun][j].c_str() ,
2681  save.wlLineList[ipPun][j] ,
2682  &relative , &absolute , LineType ) ) <=0 )
2683  {
2684  if( !h2.lgEnabled && strncmp( save.chLineListLabel[ipPun][j].c_str() , "H2 " , 4 )==0 )
2685  {
2686  static bool lgMustPrintFirstTime = true;
2687  if( lgMustPrintFirstTime )
2688  {
2689  /* it's an H2 line and H2 is not being done - ignore it */
2690  fprintf( ioQQQ,"Did not find an H2 line, the large model is not "
2691  "included, so I will ignore it. Log intensity set to -30.\n" );
2692  fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n");
2693  lgMustPrintFirstTime = false;
2694  }
2695  relative = -30.f;
2696  absolute = -30.f;
2697  }
2698  else if( lgAbort )
2699  {
2700  /* we are in abort mode */
2701  relative = -30.f;
2702  absolute = -30.f;
2703  }
2704  else
2705  {
2706  lgBadLine = true;
2707  }
2708  }
2709 
2710  /* options to do either relative or absolute intensity
2711  * default is relative, is absolute keyword on line then
2712  * punarg set to 1 */
2713  /* straight line intensities */
2714  if( save.punarg[ipPun][0] > 0 )
2715  PrtQuantity = absolute;
2716  else
2717  PrtQuantity = relative;
2718 
2719  // column mode, print label
2720  if( save.punarg[ipPun][1] )
2721  {
2722  /* if taking ratio then put div sign between pairs */
2723  if( save.lgLineListRatio[ipPun] && is_odd(j) )
2724  fprintf( save.params[ipPun].ipPnunit , "/" );
2725 
2726  fprintf( save.params[ipPun].ipPnunit, "%s ", save.chLineListLabel[ipPun][j].c_str() );
2727  char chTemp[100];
2728  sprt_wl( chTemp, save.wlLineList[ipPun][j] );
2729  fprintf( save.params[ipPun].ipPnunit, "%s ", chTemp );
2730  }
2731 
2732  /* if taking ratio print every other line as ratio
2733  * with previous line */
2734  if( save.lgLineListRatio[ipPun] )
2735  {
2736  /* do line pair ratios */
2737  static double SaveQuantity = 0;
2738  if( is_odd(j) )
2739  fprintf( save.params[ipPun].ipPnunit, "\t%.4e" ,
2740  SaveQuantity / SDIV( PrtQuantity ) );
2741  else
2742  SaveQuantity = PrtQuantity;
2743  }
2744  else
2745  {
2746  fprintf( save.params[ipPun].ipPnunit, "\t%.4e" , PrtQuantity );
2747  }
2748  // column printout, but check if first of pair
2749  if( save.punarg[ipPun][1] )
2750  {
2751  if( !save.lgLineListRatio[ipPun] ||
2752  is_odd(j) )
2753  fprintf( save.params[ipPun].ipPnunit, "\n" );
2754  }
2755  }
2756  fprintf( save.params[ipPun].ipPnunit, "\n" );
2757  if( lgBadLine )
2758  {
2759  fprintf(ioQQQ,"DISASTER - did not find line(s) in the Line List table\n");
2760  cdEXIT(EXIT_FAILURE);
2761  }
2762  }
2763  }
2764 
2765  else if( strcmp(save.chSave[ipPun],"MAP ") == 0 )
2766  {
2767  /* do the map now if we are at the zone, or if this
2768  * is the LAST call to this routine and map not done yet */
2769  if( !hcmap.lgMapDone &&
2770  (nzone == hcmap.MapZone || lgLastOnly ) )
2771  {
2772  bool lgTlkSav = called.lgTalk;
2773  called.lgTalk = cpu.i().lgMPI_talk();
2774  hcmap.lgMapBeingDone = true;
2775  map_do(save.params[ipPun].ipPnunit , " map");
2776  called.lgTalk = lgTlkSav;
2777  }
2778  }
2779 
2780  // save molecules
2781  else if( strcmp(save.chSave[ipPun],"MOLE") == 0 )
2782  {
2783  if( save.lgSaveHeader(ipPun) )
2784  {
2785  fprintf( save.params[ipPun].ipPnunit,
2786  "#depth\tAV(point)\tAV(extend)\tCO diss rate\tC recom rate");
2787 
2788  for(i=0; i<mole_global.num_calc; ++i )
2789  {
2790  if( mole_global.list[i]->n_react > 0 )
2791  fprintf( save.params[ipPun].ipPnunit, "\t%s",
2792  mole_global.list[i]->label.c_str() );
2793  }
2794  fprintf ( save.params[ipPun].ipPnunit, "\n");
2795  save.SaveHeaderDone(ipPun);
2796  }
2797  if( ! lgLastOnly )
2798  {
2799  /* molecules, especially for PDR, first give radius */
2800  fprintf( save.params[ipPun].ipPnunit, "%.5e\t" , radius.depth_mid_zone );
2801 
2802  /* visual extinction of point source (star)*/
2803  fprintf( save.params[ipPun].ipPnunit, "%.5e\t" , rfield.extin_mag_V_point);
2804 
2805  /* visual extinction of an extended source (like a PDR)*/
2806  fprintf( save.params[ipPun].ipPnunit, "%.5e\t" , rfield.extin_mag_V_extended);
2807 
2808  /* carbon monoxide photodissociation rate */
2809  fprintf( save.params[ipPun].ipPnunit, "%.5e\t" , mole.findrk("PHOTON,CO=>C,O") );
2810 
2811  /* carbon recombination rate */
2812  fprintf( save.params[ipPun].ipPnunit, "%.5e" , ionbal.RateRecomTot[ipCARBON][0] );
2813 
2814  /* now do all the molecules */
2815  for(j=0; j<mole_global.num_calc; ++j )
2816  {
2817  if( mole_global.list[j]->n_react > 0 )
2818  fprintf( save.params[ipPun].ipPnunit, "\t%.5e",
2819  mole.species[j].den );
2820  }
2821 
2822  fprintf(save.params[ipPun].ipPnunit,"\n");
2823  }
2824  }
2825 
2826  else if( strcmp(save.chSave[ipPun],"OPAC") == 0 )
2827  {
2828  /* save opacity- routine will parse which type of opacity save to do */
2829  if( save.lgSaveEveryZone[ipPun] || lgLastOnly )
2830  save_opacity(save.params[ipPun].ipPnunit,ipPun);
2831  }
2832 
2833  /* save coarse optical depths command */
2834  else if( strcmp(save.chSave[ipPun],"OPTc") == 0 )
2835  {
2836  if( save.lgSaveEveryZone[ipPun] || lgLastOnly )
2837  {
2838  for( j=0; j < rfield.nflux; j++ )
2839  {
2840  fprintf( save.params[ipPun].ipPnunit,
2841  "%13.5e\t%.3e\t%12.4e\t%.3e\n",
2842  AnuUnit(rfield.anu(j)),
2843  opac.TauAbsFace[j]+opac.TauScatFace[j],
2844  opac.TauAbsFace[j],
2845  opac.TauScatFace[j] );
2846  }
2847  }
2848  }
2849 
2850  /* save fine optical depths command */
2851  else if( strcmp(save.chSave[ipPun],"OPTf") == 0 )
2852  {
2853  if( save.lgSaveEveryZone[ipPun] || lgLastOnly )
2854  {
2855  long nu_hi , nskip;
2856  if( save.punarg[ipPun][0] > 0. )
2857  /* possible lower bounds to energy range - will be zero if not set */
2858  j = ipFineCont( save.punarg[ipPun][0] );
2859  else
2860  j = 0;
2861 
2862  /* upper limit */
2863  if( save.punarg[ipPun][1]> 0. )
2864  nu_hi = ipFineCont( save.punarg[ipPun][1]);
2865  else
2866  nu_hi = rfield.nfine;
2867 
2868  /* we will bring nskip cells together into one printed
2869  * number to make output smaller - default is 10 */
2870  nskip = (long)abs(save.punarg[ipPun][2]);
2871  nskip = MAX2( 1, nskip );
2872 
2873  do
2874  {
2875  realnum sum1 = rfield.fine_opt_depth[j];
2876  realnum sum2 = rfield.fine_opac_zone[j];
2877  /* want to report the central wavelength of the cell */
2878  realnum xnu = rfield.fine_anu[j];
2879  for( long jj=1; jj<nskip; ++jj )
2880  {
2881  sum1 += rfield.fine_opt_depth[j+jj];
2882  sum2 += rfield.fine_opac_zone[j+jj];
2883  xnu += rfield.fine_anu[j+jj];
2884  }
2885  // report each point, even 0, if ALL keyword appears
2886  if( sum2>0. || save.punarg[ipPun][2]<0)
2887  fprintf( save.params[ipPun].ipPnunit,
2888  "%12.6e\t%.3e\t%.3e\n",
2889  AnuUnit(xnu/nskip),
2890  sum1/nskip ,
2891  sum2/nskip);
2892  j += nskip;
2893  }while( j < nu_hi );
2894  }
2895  }
2896 
2897  else if( strcmp(save.chSave[ipPun]," OTS") == 0 )
2898  {
2899  double ConMax = 0.;
2900  double xLinMax = 0.;
2901  double opConSum = 0.;
2902  double opLinSum = 0.;
2903  long ipLinMax = 1;
2904  long ipConMax = 1;
2905 
2906  for( j=0; j < rfield.nflux; j++ )
2907  {
2908  opConSum += rfield.otscon[j]*opac.opacity_abs[j];
2909  opLinSum += rfield.otslin[j]*opac.opacity_abs[j];
2910  if( rfield.otslin[j]*opac.opacity_abs[j] > xLinMax )
2911  {
2912  xLinMax = rfield.otslin[j]*opac.opacity_abs[j];
2913  ipLinMax = j+1;
2914  }
2915  if( rfield.otscon[j]*opac.opacity_abs[j] > ConMax )
2916  {
2917  ConMax = rfield.otscon[j]*opac.opacity_abs[j];
2918  ipConMax = j+1;
2919  }
2920  }
2921  fprintf( save.params[ipPun].ipPnunit,
2922  "tot con lin=%.2e%.2e lin=%.4s%.4e%.2e con=%.4s%.4e%.2e\n",
2923  opConSum, opLinSum, rfield.chLineLabel[ipLinMax-1].c_str()
2924  , rfield.anu(ipLinMax-1), xLinMax, rfield.chContLabel[ipConMax-1].c_str()
2925  , rfield.anu(ipConMax-1), ConMax );
2926  }
2927 
2928  else if( strcmp(save.chSave[ipPun],"OVER") == 0 )
2929  {
2930  /* save overview
2931  * this is the floor for the smallest ionization fractions printed */
2932  double toosmall = SMALLFLOAT ,
2933  hold;
2934 
2935  /* overview of model results,
2936  * depth, te, hden, eden, ion fractions H, He, c, O */
2937  if( ! lgLastOnly )
2938  {
2939 
2940  /* print the depth */
2941  fprintf( save.params[ipPun].ipPnunit, "%.5e\t", radius.depth_mid_zone );
2942 
2943  /* temperature, heating */
2944  if(dynamics.Cool() > dynamics.Heat())
2945  {
2946  fprintf( save.params[ipPun].ipPnunit, "%.4e\t%.3e",
2947  PrtLogLin(phycon.te ),
2948  PrtLogLin(thermal.htot-dynamics.Heat() ) );
2949  }
2950  else
2951  {
2952  double diff = fabs(thermal.htot-dynamics.Cool());
2953  fprintf( save.params[ipPun].ipPnunit, "%.4e\t%.3e",
2954  PrtLogLin(phycon.te),
2955  PrtLogLin( diff ) );
2956  }
2957 
2958  /* hydrogen and electron densities */
2959  fprintf( save.params[ipPun].ipPnunit, "\t%.4e\t%.4e",
2960  PrtLogLin(dense.gas_phase[ipHYDROGEN]),
2961  PrtLogLin(dense.eden ) );
2962 
2963  /* molecular fraction of hydrogen */
2964  fprintf( save.params[ipPun].ipPnunit, "\t%.4e",
2965  PrtLogLin(MAX2(toosmall,2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN])) );
2966 
2967  /* ionization fractions of hydrogen */
2968  fprintf( save.params[ipPun].ipPnunit, "\t%.4e\t%.4e",
2969  PrtLogLin(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])),
2970  PrtLogLin(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])) );
2971 
2972  /* ionization fractions of helium */
2973  for( j=1; j <= 3; j++ )
2974  {
2975  double arg1 = SDIV(dense.gas_phase[ipHELIUM]);
2976  arg1 = MAX2(toosmall,dense.xIonDense[ipHELIUM][j-1]/arg1 );
2977  fprintf( save.params[ipPun].ipPnunit, "\t%.4e",
2978  PrtLogLin(arg1) );
2979  }
2980 
2981  /* carbon monoxide molecular fraction of CO */
2982  hold = SDIV(dense.gas_phase[ipCARBON]);
2983  hold = findspecieslocal("CO")->den/hold;
2984  hold = MAX2(toosmall, hold );
2985  fprintf( save.params[ipPun].ipPnunit, "\t%.4e", PrtLogLin(hold) );
2986 
2987  /* ionization fractions of carbon */
2988  for( j=1; j <= 4; j++ )
2989  {
2990  hold = SDIV(dense.gas_phase[ipCARBON]);
2991  hold = MAX2(toosmall,dense.xIonDense[ipCARBON][j-1]/hold);
2992  fprintf( save.params[ipPun].ipPnunit, "\t%.4e",
2993  PrtLogLin(hold) );
2994  }
2995 
2996  /* ionization fractions of oxygen */
2997  for( j=1; j <= 6; j++ )
2998  {
2999  hold = SDIV(dense.gas_phase[ipOXYGEN]);
3000  hold = MAX2(toosmall,dense.xIonDense[ipOXYGEN][j-1]/hold);
3001  fprintf( save.params[ipPun].ipPnunit, "\t%.4e",
3002  PrtLogLin(hold) );
3003  }
3004 
3005  // molecular fraction of H2O
3006  hold = SDIV(dense.gas_phase[ipOXYGEN]);
3007  hold = findspecieslocal("H2O")->den/hold;
3008  hold = MAX2(toosmall, hold );
3009  fprintf( save.params[ipPun].ipPnunit, "\t%.4e", PrtLogLin(hold) );
3010 
3011  /* visual extinction of point source (star)*/
3012  fprintf( save.params[ipPun].ipPnunit, "\t%.2e" , rfield.extin_mag_V_point);
3013 
3014  /* visual extinction of an extended source (like a PDR)*/
3015  fprintf( save.params[ipPun].ipPnunit, "\t%.2e" , rfield.extin_mag_V_extended);
3016 
3017  /* >>chng 20 07 25 add Opt Depth @ 1 Ryd as per Vianney LeBouteiller post */
3018  /*Lyman continuum optical depth to illuminated face */
3019  fprintf( save.params[ipPun].ipPnunit, "\t%.4e\n", opac.TauAbsFace[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon]);
3020  }
3021  }
3022 
3023  else if( strcmp(save.chSave[ipPun]," PDR") == 0 )
3024  {
3025  /* this is the save PDR command */
3026  if( ! lgLastOnly )
3027  {
3028  /* convert optical depth at wavelength of V filter
3029  * into magnitudes of extinction */
3030  /* >>chyng 03 feb 25, report extinction to illuminated face,
3031  * rather than total extinction which included far side when
3032  * sphere was set */
3033  /*av = opac.TauTotalGeo[0][rfield.ipV_filter-1]*1.08574;*/
3034 
3035  fprintf( save.params[ipPun].ipPnunit,
3036  "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t",
3037  radius.depth_mid_zone,
3038  /* total hydrogen column density, all forms */
3039  colden.colden[ipCOL_HTOT],
3040  phycon.te,
3041  /* fraction of H that is atomic */
3042  dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
3043  /* ratio of n(H2) to total H, == 0.5 when fully molecular */
3044  2.*findspecieslocal("H2")->den/dense.gas_phase[ipHYDROGEN],
3045  2.*findspecieslocal("H2*")->den/dense.gas_phase[ipHYDROGEN],
3046  /* atomic to total carbon */
3047  dense.xIonDense[ipCARBON][0]/SDIV(dense.gas_phase[ipCARBON]),
3048  findspecieslocal("CO")->den/SDIV(dense.gas_phase[ipCARBON]),
3049  findspecieslocal("H2O")->den/SDIV(dense.gas_phase[ipOXYGEN]),
3050  /* hmi.UV_Cont_rel2_Habing_TH85 is field relative to Habing background, dimensionless */
3051  hmi.UV_Cont_rel2_Habing_TH85_depth);
3052 
3053  /* visual extinction due to dust alone, of point source (star)*/
3054  fprintf( save.params[ipPun].ipPnunit, "%.2e\t" , rfield.extin_mag_V_point);
3055 
3056  /* visual extinction due to dust alone, of an extended source (like a PDR)*/
3057  fprintf( save.params[ipPun].ipPnunit, "%.2e\t" , rfield.extin_mag_V_extended);
3058 
3059  /* visual extinction (all sources) of a point source (like a PDR)*/
3060  fprintf( save.params[ipPun].ipPnunit, "%.2e\n", opac.TauAbsGeo[0][rfield.ipV_filter] );
3061  }
3062  }
3063 
3064  /* performance characteristics per zone */
3065  else if( strcmp(save.chSave[ipPun],"PERF") == 0 )
3066  {
3067  if( save.lgSaveHeader(ipPun) )
3068  {
3069  fprintf( save.params[ipPun].ipPnunit,
3070  "#zone\tdTime\tElapsed t\tnPres2Ioniz" );
3071  for( size_t i = 0; i < conv.ntypes(); i++ )
3072  {
3073  fprintf( save.params[ipPun].ipPnunit, "\t%s",
3074  conv.getCounterName(i) );
3075  }
3076  fprintf( save.params[ipPun].ipPnunit, "\n" );
3077  save.SaveHeaderDone(ipPun);
3078  }
3079  if( ! lgLastOnly )
3080  {
3081  static double ElapsedTime , ZoneTime;
3082  if( nzone<=1 )
3083  {
3084  ElapsedTime = cdExecTime();
3085  ZoneTime = 0.;
3086  }
3087  else
3088  {
3089  double t = cdExecTime();
3090  ZoneTime = t - ElapsedTime;
3091  ElapsedTime = t;
3092  }
3093 
3094  /* zone, time for this zone, elapsed time */
3095  fprintf( save.params[ipPun].ipPnunit, " %ld\t%.3f\t%.2f\t%li",
3096  nzone, ZoneTime , ElapsedTime, conv.nPres2Ioniz );
3097  // print various loop counters
3098  for( size_t i=0; i<conv.ntypes(); ++i )
3099  fprintf( save.params[ipPun].ipPnunit, "\t%li", conv.getCounterZone(i) );
3100  fprintf( save.params[ipPun].ipPnunit, "\n" );
3101  }
3102  }
3103 
3104  else if( strcmp(save.chSave[ipPun],"PHYS") == 0 )
3105  {
3106  if( ! lgLastOnly )
3107  {
3108  /* save physical conditions */
3109  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.4e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
3110  radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN],
3111  dense.eden, thermal.htot, wind.AccelTotalOutward, geometry.FillFac );
3112  }
3113  }
3114 
3115  else if( strcmp(save.chSave[ipPun],"PRES") == 0 )
3116  {
3117  /* the save pressure command */
3118  if( ! lgLastOnly )
3119  {
3120  fprintf( save.params[ipPun].ipPnunit,
3121  "%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\t%c\n",
3122  /*A 1 #P depth */
3123  radius.depth_mid_zone,
3124  /*B 2 Perror */
3125  pressure.PresTotlError*100.,
3126  /*C 3 Pcurrent */
3127  pressure.PresTotlCurr,
3128  /*D 4 Pln + pintg
3129  * >>chng 06 apr 19, subtract pinzon the acceleration added in this zone
3130  * since is not total at outer edge of zone, above is at inner edge */
3131  pressure.PresTotlInit + pressure.PresInteg - pressure.pinzon,
3132  /*E 5 pgas (0) */
3133  pressure.PresTotlInit,
3134  /*F 6 Pgas */
3135  pressure.PresGasCurr,
3136  /*G 7 Pram */
3137  pressure.PresRamCurr,
3138  /*H 8 P rad in lines */
3139  pressure.pres_radiation_lines_curr,
3140  /*I 9 Pinteg subtract continuum rad pres which has already been added on */
3141  pressure.PresInteg - pressure.pinzon,
3142  /*J 10 V(wind km/s) wind speed in km/s */
3143  wind.windv/1e5,
3144  /*K cad(km/s) sound speed in km/s */
3145  timesc.sound_speed_adiabatic/1e5,
3146  /* the magnetic pressure */
3147  magnetic.pressure ,
3148  /* the local turbulent velocity in km/s */
3149  DoppVel.TurbVel/1e5 ,
3150  /* turbulent pressure */
3151  pressure.PresTurbCurr*DoppVel.lgTurb_pressure ,
3152  /* gravitational pressure */
3153  pressure.IntegRhoGravity,
3154  // the integral of electron scattering acceleration in
3155  // the absence of any absorptio, minus acceleration in current
3156  // zone, which has been added in - done this way, result is
3157  // zero in first zonen
3158  pressure.PresIntegElecThin-pressure.pinzon_PresIntegElecThin,
3159  // is this converged?
3160  TorF(conv.lgConvPres) );
3161  }
3162  }
3163  else if( strcmp(save.chSave[ipPun],"PREL") == 0 )
3164  {
3165  /* line pressure contributors */
3166  fprintf( save.params[ipPun].ipPnunit,
3167  "%.5e\t%.3e\t%.3e\t",
3168  /*A 1 #P depth */
3169  radius.depth_mid_zone ,
3170  pressure.PresTotlCurr,
3171  pressure.pres_radiation_lines_curr/SDIV(pressure.PresTotlCurr) );
3172  PrtLinePres(save.params[ipPun].ipPnunit);
3173 
3174  }
3175 
3176  else if( save.chSave[ipPun][0]=='R' )
3177  {
3178  /* work around internal limits to Microsoft vs compiler */
3179  if( strcmp(save.chSave[ipPun],"RADI") == 0 )
3180  {
3181  /* save radius information for all zones */
3182  if( ! lgLastOnly )
3183  {
3184  fprintf( save.params[ipPun].ipPnunit, "%ld\t%.5e\t%.4e\t%.4e\n",
3185  nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
3186  radius.drad );
3187  }
3188  }
3189 
3190  else if( strcmp(save.chSave[ipPun],"RADO") == 0 )
3191  {
3192  /* save radius information for only the last zone */
3193  if( lgLastOnly )
3194  {
3195  fprintf( save.params[ipPun].ipPnunit, "%ld\t%.5e\t%.4e\t%.4e\n",
3196  nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
3197  radius.drad );
3198  }
3199  }
3200 
3201  else if( strcmp(save.chSave[ipPun],"RESU") == 0 )
3202  {
3203  /* save results of the calculation */
3204  if( lgLastOnly )
3205  SaveResults(save.params[ipPun].ipPnunit);
3206  }
3207 
3208  else if( strcmp(save.chSave[ipPun],"RECA") == 0 )
3209  {
3210  /* this will create table for AGN3 then exit,
3211  * routine is in makerecom.c */
3212  ion_recombAGN( save.params[ipPun].ipPnunit );
3213  cdEXIT(EXIT_SUCCESS);
3214  }
3215 
3216  else if( strcmp(save.chSave[ipPun],"RECE") == 0 )
3217  {
3218  /* save recombination efficiencies,
3219  * option turned on with the "save recombination efficiencies" command
3220  * output for the save recombination coefficients command is actually
3221  * produced by a series of routines, as they generate the recombination
3222  * coefficients. these include
3223  * dielsupres, helium, hydrorecom, iibod, and makerecomb*/
3224  fprintf( save.params[ipPun].ipPnunit,
3225  "%12.4e %12.4e %12.4e %12.4e\n",
3226  iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].RadRecomb[ipRecRad],
3227  iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].RadRecomb[ipRecNetEsc] ,
3228  iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].RadRecomb[ipRecRad],
3229  iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].RadRecomb[ipRecNetEsc]);
3230  }
3231 
3232  else
3233  {
3234  /* this can't happen */
3235  TotalInsanity();
3236  }
3237  }
3238 
3239  else if( strcmp(save.chSave[ipPun],"SECO") == 0 )
3240  {
3241  /* save secondary ionization */
3242  if( ! lgLastOnly )
3243  fprintf(save.params[ipPun].ipPnunit,
3244  "%.5e\t%.3e\t%.3e\t%.3e\n",
3245  radius.depth ,
3246  secondaries.csupra[ipHYDROGEN][0],
3247  secondaries.csupra[ipHYDROGEN][0]*2.02,
3248  secondaries.x12tot );
3249  }
3250 
3251  else if( strcmp(save.chSave[ipPun],"SOUS") == 0 )
3252  {
3253  /* full spectrum of continuum source function at 1 depth
3254  * command was "save source spectrum" */
3255  if( ! lgLastOnly )
3256  {
3257  long limit = MIN2(rfield.ipMaxBolt,rfield.nflux);
3258  for( j=0; j < limit; j++ )
3259  {
3260  fprintf( save.params[ipPun].ipPnunit,
3261  "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
3262  AnuUnit(rfield.anu(j)),
3263  rfield.ConEmitLocal[nzone][j]/rfield.widflx(j),
3264  opac.opacity_abs[j],
3265  rfield.ConSourceFcnLocal[nzone][j],
3266  rfield.ConSourceFcnLocal[nzone][j]/plankf(j),
3267  safe_div(rfield.ConSourceFcnLocal[nzone][j],rfield.flux[0][j]));
3268  }
3269  }
3270  }
3271 
3272  else if( strcmp(save.chSave[ipPun],"SOUD") == 0 )
3273  {
3274  /* parts of continuum source function vs depth
3275  * command was save source function depth */
3276  j = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon + 2;
3277  fprintf( save.params[ipPun].ipPnunit,
3278  "%.4e\t%.4e\t%.4e\t%.4e\n",
3279  opac.TauAbsFace[j-1],
3280  rfield.ConEmitLocal[nzone][j-1]/rfield.widflx(j-1)/MAX2(1e-35,opac.opacity_abs[j-1]),
3281  rfield.otscon[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1],
3282  rfield.otscon[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/opac.opacity_abs[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] );
3283  }
3284 
3285  /* this is save special option */
3286  else if( strcmp(save.chSave[ipPun],"SPEC") == 0 )
3287  {
3288  SaveSpecial(save.params[ipPun].ipPnunit,chTime);
3289  }
3290 
3291  /* this is save species option */
3292  else if( strcmp(save.chSave[ipPun],"SPCS") == 0 )
3293  {
3294  if( strncmp( save.chSaveArgs[ipPun], "CON", 3 ) == 0 )
3295  {
3296  if( lgLastOnly )
3297  SaveSpeciesPseudoCont( ipPun,
3298  save.chSaveSpecies[ipPun][0] );
3299  }
3300  else if( strcmp( save.chSaveArgs[ipPun], "BAND" ) == 0 )
3301  {
3302  if( lgLastOnly )
3303  SaveSpeciesBands( ipPun,
3304  save.chSaveSpecies[ipPun][0],
3305  save.SpeciesBandFile[ipPun] );
3306  }
3307  else if( strcmp( save.chSaveArgs[ipPun], "OPTD" ) == 0 )
3308  {
3309  if( lgLastOnly )
3310  SaveSpeciesOptDep( ipPun, save.chSaveSpecies[ipPun][0] );
3311  }
3312  else if( ( ! lgLastOnly && strcmp(save.chSaveArgs[ipPun],"COLU") != 0 ) ||
3313  ( lgLastOnly && strcmp(save.chSaveArgs[ipPun],"COLU") == 0 ) )
3314  SaveSpecies(save.params[ipPun].ipPnunit , ipPun);
3315  }
3316 
3317  else if( strcmp(save.chSave[ipPun],"TEMP") == 0 )
3318  {
3319  static double deriv_old=-1;
3320  double deriv=-1. , deriv_sec;
3321  /* temperature and its derivatives */
3322  fprintf( save.params[ipPun].ipPnunit, "%.5e\t%.4e\t%.2e",
3323  radius.depth_mid_zone,
3324  phycon.te,
3325  thermal.dCooldT );
3326  /* if second zone then have one deriv */
3327  if( nzone >1 )
3328  {
3329  deriv = (phycon.te - struc.testr[nzone-2])/ radius.drad;
3330  fprintf( save.params[ipPun].ipPnunit, "\t%.2e", deriv );
3331  /* if third zone then have second deriv */
3332  if( nzone > 2 )
3333  {
3334  deriv_sec = (deriv-deriv_old)/ radius.drad;
3335  fprintf( save.params[ipPun].ipPnunit, "\t%.2e",
3336  deriv_sec );
3337  }
3338  deriv_old = deriv;
3339  }
3340  fprintf( save.params[ipPun].ipPnunit, "\n");
3341  }
3342 
3343  /* time dependent model */
3344  else if( strcmp(save.chSave[ipPun],"TIMD") == 0 )
3345  {
3346  if( lgLastOnly )
3347  DynaPunchTimeDep( save.params[ipPun].ipPnunit , "END" );
3348  }
3349 
3350  /* execution time per zone */
3351  else if( strcmp(save.chSave[ipPun],"XTIM") == 0 )
3352  {
3353  static double ElapsedTime , ZoneTime;
3354  if( nzone<=1 )
3355  {
3356  ElapsedTime = cdExecTime();
3357  ZoneTime = 0.;
3358  }
3359  else
3360  {
3361  double t = cdExecTime();
3362  ZoneTime = t - ElapsedTime;
3363  ElapsedTime = t;
3364  }
3365 
3366  /* zone, time for this zone, elapsed time */
3367  fprintf( save.params[ipPun].ipPnunit, " %ld\t%.3f\t%.2f\n",
3368  nzone, ZoneTime , ElapsedTime );
3369  }
3370 
3371  else if( strcmp(save.chSave[ipPun],"TPRE") == 0 )
3372  {
3373  /* temperature and its predictors, turned on with save tprid */
3374  fprintf( save.params[ipPun].ipPnunit, "%5ld %11.4e %11.4e %11.4e %g\n",
3375  nzone, phycon.TeInit, phycon.TeProp, phycon.te,
3376  (phycon.TeProp- phycon.te)/phycon.te );
3377  }
3378 
3379  else if( strcmp(save.chSave[ipPun],"WIND") == 0 )
3380  {
3381  /* wind velocity, radiative acceleration, and ratio total
3382  * to electron scattering acceleration */
3383  /* first test only save last zone */
3384  if( (save.punarg[ipPun][0] == 0 && lgLastOnly)
3385  ||
3386  /* this test save all zones */
3387  (save.punarg[ipPun][0] == 1 && ! lgLastOnly ) )
3388  {
3389  fprintf( save.params[ipPun].ipPnunit,
3390  "%.5e\t%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
3391  radius.Radius_mid_zone,
3392  radius.depth_mid_zone,
3393  wind.windv,
3394  wind.AccelTotalOutward,
3395  wind.AccelLine,
3396  wind.AccelCont ,
3397  wind.fmul ,
3398  wind.AccelGravity );
3399  }
3400  }
3401 
3402  else if( strcmp(save.chSave[ipPun],"XATT") == 0 )
3403  {
3404  /* attenuated incident continuum */
3405  ASSERT( grid.lgOutputTypeOn[2] );
3406 
3407  if( lgLastOnly )
3408  saveFITSfile( save.params[ipPun].ipPnunit, 2 );
3409  }
3410  else if( strcmp(save.chSave[ipPun],"XRFI") == 0 )
3411  {
3412  /* reflected incident continuum */
3413  ASSERT( grid.lgOutputTypeOn[3] );
3414 
3415  if( lgLastOnly )
3416  saveFITSfile( save.params[ipPun].ipPnunit, 3 );
3417  }
3418  else if( strcmp(save.chSave[ipPun],"XINC") == 0 )
3419  {
3420  /* incident continuum */
3421  ASSERT( grid.lgOutputTypeOn[1] );
3422 
3423  if( lgLastOnly )
3424  saveFITSfile( save.params[ipPun].ipPnunit, 1 );
3425  }
3426  else if( strcmp(save.chSave[ipPun],"XDFR") == 0 )
3427  {
3428  /* reflected diffuse continuous emission */
3429  ASSERT( grid.lgOutputTypeOn[5] );
3430 
3431  if( lgLastOnly )
3432  saveFITSfile( save.params[ipPun].ipPnunit, 5 );
3433  }
3434  else if( strcmp(save.chSave[ipPun],"XDFO") == 0 )
3435  {
3436  /* diffuse continuous emission outward */
3437  ASSERT( grid.lgOutputTypeOn[4] );
3438 
3439  if( lgLastOnly )
3440  saveFITSfile( save.params[ipPun].ipPnunit, 4 );
3441  }
3442  else if( strcmp(save.chSave[ipPun],"XLNR") == 0 )
3443  {
3444  /* reflected lines */
3445  ASSERT( grid.lgOutputTypeOn[7] );
3446 
3447  if( lgLastOnly )
3448  saveFITSfile( save.params[ipPun].ipPnunit, 7 );
3449  }
3450  else if( strcmp(save.chSave[ipPun],"XLNO") == 0 )
3451  {
3452  /* outward lines */
3453  ASSERT( grid.lgOutputTypeOn[6] );
3454 
3455  if( lgLastOnly )
3456  saveFITSfile( save.params[ipPun].ipPnunit, 6 );
3457  }
3458  else if( strcmp(save.chSave[ipPun],"XREF") == 0 )
3459  {
3460  /* total reflected, lines and continuum */
3461  ASSERT( grid.lgOutputTypeOn[9] );
3462 
3463  if( lgLastOnly )
3464  saveFITSfile( save.params[ipPun].ipPnunit, 9 );
3465  }
3466  else if( strcmp(save.chSave[ipPun],"XTOT") == 0 )
3467  {
3468  /* total spectrum, reflected plus transmitted */
3469  ASSERT( grid.lgOutputTypeOn[0] );
3470 
3471  if( lgLastOnly )
3472  saveFITSfile( save.params[ipPun].ipPnunit, 0 );
3473  }
3474  else if( strcmp(save.chSave[ipPun],"XTRN") == 0 )
3475  {
3476  /* total outward, lines and continuum */
3477  ASSERT( grid.lgOutputTypeOn[8] );
3478 
3479  if( lgLastOnly )
3480  saveFITSfile( save.params[ipPun].ipPnunit, 8 );
3481  }
3482  else if( strcmp(save.chSave[ipPun],"XSPM") == 0 )
3483  {
3484  /* exp(-tau) to the illuminated face */
3485  ASSERT( grid.lgOutputTypeOn[10] );
3486 
3487  if( lgLastOnly )
3488  saveFITSfile( save.params[ipPun].ipPnunit, 10 );
3489  }
3490  // termination of second set of nested if's
3491  // error if we have not matched key
3492  /* there are a few "save" commands that are handled elsewhere
3493  * save dr is an example. These will have lgRealSave set false */
3494  // lgNoHitFirstBranch says did not find in previous nest of if's
3495  else if( save.lgRealSave[ipPun] && lgNoHitFirstBranch )
3496  {
3497  /* this is insanity, internal flag set in ParseSave not seen here */
3498  fprintf( ioQQQ, " PROBLEM DISASTER SaveDo does not recognize flag %4.4s set by ParseSave. This is impossible.\n",
3499  save.chSave[ipPun] );
3500  TotalInsanity();
3501  }
3502 
3503  /* print special hash string to separate out various iterations
3504  * chTime is LAST on last iteration
3505  * save.lgHashEndIter flag is true by default, set false
3506  * with "no hash" keyword on save command
3507  * save.lg_separate_iterations is true by default, set false
3508  * when save time dependent calcs since do not want special
3509  * character between time steps
3510  * grid.lgGrid is only true when doing a grid of calculations */
3511  if( lgLastOnly &&
3512  !(iterations.lgLastIt && !grid.lgGrid ) &&
3513  save.lgHashEndIter[ipPun] &&
3514  save.lg_separate_iterations[ipPun] &&
3515  !save.lgFITS[ipPun] )
3516  {
3517  if( dynamics.lgTimeDependentStatic && strcmp( save.chHashString , "TIME_DEP" )==0 )
3518  {
3519  fprintf( save.params[ipPun].ipPnunit, "\"time=%f\n",
3520  dynamics.time_elapsed );
3521  }
3522  else if( strcmp( save.chHashString , "\n" )==0 )
3523  {
3524  fprintf( save.params[ipPun].ipPnunit, "%s\n",
3525  save.chHashString );
3526  }
3527  else
3528  {
3529  fprintf( save.params[ipPun].ipPnunit, "%s",
3530  save.chHashString );
3531  if( grid.lgGrid && ( iterations.lgLastIt || lgAbort ) )
3532  fprintf( save.params[ipPun].ipPnunit, " GRID_DELIMIT -- grid%09ld",
3533  optimize.nOptimiz );
3534  fprintf( save.params[ipPun].ipPnunit, "\n" );
3535  }
3536  }
3537  if( save.lgFLUSH )
3538  fflush( save.params[ipPun].ipPnunit );
3539  }
3540  }
3541  return;
3542 }
3543 
3544 /*SaveLineIntensity produce the 'save lines intensity' output */
3545 STATIC void SaveLineIntensity(FILE * ioPUN, long int ipPun , realnum Threshold )
3546 {
3547  long int i;
3548 
3549  DEBUG_ENTRY( "SaveLineIntensity()" );
3550 
3551  /* used to save out all the emission line intensities
3552  * first initialize the line image reader */
3553 
3554  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3555  input.echo(ioPUN);
3556 
3557  /* now print any cautions or warnings */
3558  cdWarnings( ioPUN);
3559  cdCautions( ioPUN);
3560  fprintf( ioPUN, "zone=%5ld\n", nzone );
3561  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3562  fprintf( ioPUN, "begin emission lines\n" );
3563 
3564 
3565  // check whether intrinsic or emergent line emissivity
3566  bool lgEmergent = false;
3567  if( save.punarg[ipPun][0] > 0 )
3568  lgEmergent = true;
3569 
3570  /* only save non-zero intensities */
3571  fixit("value of lgEmergent isn't consistent with lgEmergent");
3572  SaveLineResults slr(ioPUN,save.lgEmergent[ipPun]);
3573 
3574  for( i=0; i < LineSave.nsum; i++ )
3575  {
3576  // Threshold is zero by default on save line intensity,
3577  // all option sets to negative number so that we report all lines
3578  if( LineSave.lines[i].SumLine(lgEmergent) > Threshold )
3579  {
3580  slr.save(&LineSave.lines[i]);
3581  }
3582  }
3583 
3584  slr.flush();
3585 
3586  fprintf( ioPUN, " \n" );
3587  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3588 
3589  return;
3590 }
3591 
3592 /* lgSaveOpticalDepths true says save optical depths */
3593 static bool lgPopsFirstCall , lgSaveOpticalDepths;
3594 
3595 /*SaveLineStuff save optical depths or source functions for all transferred lines */
3596 STATIC void SaveLineStuff(
3597  FILE * ioPUN,
3598  const char *chJob ,
3599  realnum xLimit )
3600 {
3601  DEBUG_ENTRY( "SaveLineStuff()" );
3602 
3603  /*find out which job this is and set a flag to use later */
3604  if( strcmp( &*chJob , "optical" ) == 0 )
3605  {
3606  /* save line optical depths */
3607  lgSaveOpticalDepths = true;
3608  lgPopsFirstCall = false;
3609  }
3610  else if( strcmp( &*chJob , "populat" ) == 0 )
3611  {
3612  static bool lgFirst=true;
3613  lgSaveOpticalDepths = false;
3614  /* level population information */
3615  if( lgFirst )
3616  {
3617  lgPopsFirstCall = true;
3618  fprintf(ioPUN,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n");
3619  lgFirst = false;
3620  }
3621  else
3622  {
3623  lgPopsFirstCall = false;
3624  }
3625  }
3626  else
3627  {
3628  fprintf( ioQQQ, " insane job in SaveLineStuff =%s\n",
3629  &*chJob );
3630  cdEXIT(EXIT_FAILURE);
3631  }
3632 
3633  long index = 0;
3634  /* loop over all lines, calling put1Line to create info (routine located below) */
3635  /* hydrogen like lines */
3636  /* >>chng 02 may 16, had been explicit H and He-like loops */
3637  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
3638  {
3639  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
3640  {
3641  if( dense.lgElmtOn[nelem] )
3642  {
3643  /* 06 aug 28, from numLevels_max to _local. */
3644  for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
3645  {
3646  for( long ipLo=0; ipLo <ipHi; ipLo++ )
3647  {
3648  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
3649  continue;
3650 
3651  ++index;
3652  Save1Line( iso_sp[ipISO][nelem].trans(ipHi,ipLo), ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) );
3653  }
3654  }
3655  /* also do extra Lyman lines if optical depths are to be done,
3656  * these are line that are included only for absorption, not in the
3657  * model atoms */
3658  if( lgSaveOpticalDepths )
3659  {
3660  /* >>chng 02 aug 23, for he-like, had starting on much too high a level since
3661  * index was number of levels - caught by Adrian Turner */
3662  /* now output extra line lines, starting one above those already done above */
3663  /*for( ipHi=iso_sp[ipISO][nelem].numLevels_max; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )*/
3664  /* 06 aug 28, from numLevels_max to _local. */
3665  for( long ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
3666  {
3667  ++index;
3668  Save1Line( ExtraLymanLines[ipISO][nelem][ipExtraLymanLines[ipISO][nelem][ipHi]], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) );
3669  }
3670  }
3671  }
3672  }
3673  }
3674 
3675  for( long i=0; i < nWindLine; i++ )
3676  {
3677  if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
3678  {
3679  ++index;
3680  Save1Line( TauLine2[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
3681  }
3682  }
3683 
3684  for( size_t i=0; i < UTALines.size(); i++ )
3685  {
3686  ++index;
3687  Save1Line( UTALines[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) );
3688  }
3689 
3690  /* save optical depths of H2 lines */
3691  h2.H2_PunchLineStuff( ioPUN , xLimit , index);
3692 
3693  /* data base lines */
3694  for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
3695  {
3696  realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
3697  for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
3698  em != dBaseTrans[ipSpecies].Emis().end(); ++em)
3699  {
3700  ++index;
3701  Save1Line( (*em).Tran(), ioPUN, xLimit, index, DopplerWidth );
3702  }
3703  }
3704 
3705  /*fprintf(ioPUN, "##################################\n"); */
3706  fprintf( ioPUN , "%s\n",save.chHashString );
3707  return;
3708 }
3709 
3710 /*Save1Line called by SaveLineStuff to produce output for one line */
3711 void Save1Line( const TransitionProxy& t , FILE * ioPUN , realnum xLimit , long index, realnum DopplerWidth )
3712 {
3713 
3714  if( lgSaveOpticalDepths )
3715  {
3716  /* optical depths, no special first time, only print them */
3717  if( t.Emis().TauIn() >= xLimit )
3718  {
3719  /* label like "C 4" or "H 1" */
3720  fprintf( ioPUN, "%-*.*s\t",CHARS_SPECIES, CHARS_SPECIES, chIonLbl(t).c_str());
3721 
3722  /* print wavelengths, either line in main printout labels,
3723  * or in various units in exponential notation - prt_wl is in prt.c */
3724  if( strcmp( save.chConSavEnr[save.ipConPun], "labl" )== 0 )
3725  {
3726  prt_wl( ioPUN , t.WLAng() );
3727  }
3728  else
3729  {
3730  /* this converts energy in Rydbergs into any of the other units */
3731  fprintf( ioPUN , "%.7e", AnuUnit((realnum)(t.EnergyRyd())) );
3732  }
3733  /* print the optical depth */
3734  fprintf( ioPUN , "\t%.3f", t.Emis().TauIn()*SQRTPI );
3735  /* damping constant */
3736  fprintf(ioPUN, "\t%.3e",
3737  t.Emis().dampXvel() / DopplerWidth );
3738  fprintf(ioPUN, "\n");
3739  }
3740  }
3741  else if( lgPopsFirstCall )
3742  {
3743  /* first call to line populations, print atomic parameters and indices */
3744  fprintf(ioPUN, "%li\t%s" , index , chLineLbl(t).c_str());
3745  /* stat weights */
3746  fprintf(ioPUN, "\t%.0f\t%.0f",
3747  (*t.Lo()).g() ,(*t.Hi()).g());
3748  /* energy difference, gf */
3749  fprintf(ioPUN, "\t%.2f\t%.3e",
3750  t.EnergyWN() ,t.Emis().gf());
3751  fprintf(ioPUN, "\n");
3752  }
3753  else
3754  {
3755  /* not first call, so do level populations and indices defined above */
3756  if( (*t.Hi()).Pop() > xLimit )
3757  {
3758  /* >>chng 05 may 08, add abundances, which for iso-seq species is
3759  * the density of the parent ion, for other lines, is unity.
3760  * had not been included so pops for iso seq were rel to parent ion.
3761  * caught by John Everett */
3762  /* multiplication by abundance no longer necessary since iso pops now denormalized */
3763  fprintf(ioPUN,"%li\t%.2e\t%.2e\n", index, (*t.Lo()).Pop(), (*t.Hi()).Pop() );
3764  }
3765  }
3766 }
3767 
3768 /* save AGN3 hemiss, for Chapter 4, routine is below */
3769 STATIC void AGN_Hemis(FILE *ioPUN )
3770 {
3771  const int NTE = 4;
3772  realnum te[NTE] = {5000., 10000., 15000., 20000.};
3773  realnum *agn_continuum[NTE];
3774  double TempSave = phycon.te;
3775  long i , j;
3776 
3777  DEBUG_ENTRY( "AGN_Hemis()" );
3778 
3779  /* make table of continuous emission at various temperatuers */
3780  /* first allocate space */
3781  for( i=0;i<NTE; ++i)
3782  {
3783  agn_continuum[i] = (realnum *)MALLOC((unsigned)rfield.nflux*sizeof(realnum) );
3784 
3785  /* set the next temperature */
3786  /* recompute everything at this new temp */
3787  TempChange(te[i] , true);
3788  /* converge the pressure-temperature-ionization solution for this zone */
3789  ConvPresTempEdenIoniz();
3790 
3791  /* now get the thermal emission */
3792  RT_diffuse();
3793  for(j=0;j<rfield.nflux; ++j )
3794  {
3795  agn_continuum[i][j] = rfield.ConEmitLocal[nzone][j]/(realnum)dense.eden/
3796  (dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHELIUM][1] + dense.xIonDense[ipHELIUM][2] );
3797  }
3798  }
3799 
3800  /* print title for line */
3801  fprintf(ioPUN,"wl");
3802  for( i=0;i<NTE; ++i)
3803  {
3804  fprintf(ioPUN,"\tT=%.0f",te[i]);
3805  }
3806  fprintf( ioPUN , "\tcont\n");
3807 
3808  /* not print all n temperatures across a line */
3809  for(j=0;j<rfield.nflux; ++j )
3810  {
3811  fprintf( ioPUN , "%12.5e",
3812  AnuUnit(rfield.anu(j)) );
3813  /* loop over the temperatures, and for each, calculate a continuum */
3814  for( i=0;i<NTE; ++i)
3815  {
3816  fprintf(ioPUN,"\t%.3e",agn_continuum[i][j]*rfield.anu2(j)*EN1RYD/rfield.widflx(j));
3817  }
3818  /* cont label and end of line*/
3819  fprintf( ioPUN , "\t%s\n" , rfield.chContLabel[j].c_str());
3820  }
3821 
3822  /* now free the continua */
3823  for( i=0;i<NTE; ++i)
3824  {
3825  free( agn_continuum[i] );
3826  }
3827 
3828  /* Restore temperature stored before this routine was called */
3829  /* and force update */
3830  TempChange(TempSave , true);
3831 
3832  fprintf( ioQQQ, "AGN_Hemis - result of save AGN3 hemis - I have left the code in a disturbed state, and will now exit.\n");
3833  cdEXIT(EXIT_SUCCESS);
3834 }
3835 
3836 /*SaveResults save results from save results command */
3837 /*SaveResults1Line do single line of output for the save results and save line intensity commands */
3838 STATIC void SaveResults(FILE* ioPUN)
3839 {
3840  long int i , nelem , ion;
3841 
3842  DEBUG_ENTRY( "SaveResults()" );
3843 
3844  /* used to save out line intensities, optical depths,
3845  * and column densities */
3846 
3847  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3848  input.echo(ioPUN);
3849 
3850  /* first print any cautions or warnings */
3851  cdWarnings(ioPUN);
3852  cdCautions(ioPUN);
3853  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3854 
3855  fprintf( ioPUN, "C*OPTICAL DEPTHS ELECTRON=%10.3e\n", opac.telec );
3856 
3857  fprintf( ioPUN, "BEGIN EMISSION LINES\n" );
3858  SaveLineResults slr(ioPUN,0);
3859 
3860  for( i=0; i < LineSave.nsum; i++ )
3861  {
3862  if( LineSave.lines[i].SumLine(0) > 0. )
3863  {
3864  slr.save(&LineSave.lines[i]);
3865  }
3866  }
3867 
3868  slr.flush();
3869 
3870  fprintf( ioPUN, " \n" );
3871 
3872  fprintf( ioPUN, "BEGIN COLUMN DENSITIES\n" );
3873 
3874  /* this dumps out the whole array,*/
3875  /* following loop relies on LIMELM being 30, assert it here in case
3876  * this is ever changed */
3877  ASSERT( LIMELM == 30 );
3878  /* this order of indices is to keep 30 as the fastest variable,
3879  * and the 32 (LIMELM+1) as the slower one */
3880  for( nelem=0; nelem<LIMELM; nelem++ )
3881  {
3882  for(ion=0; ion < nelem+1; ion++)
3883  {
3884  fprintf( ioPUN, " %10.3e", mean.xIonMean[0][nelem][ion][0] );
3885  /* throw line feed every 10 numbers */
3886  if( nelem==9|| nelem==19 || nelem==29 )
3887  {
3888  fprintf( ioPUN, "\n" );
3889  }
3890  }
3891  }
3892 
3893  fprintf( ioPUN, "END OF RESULTS\n" );
3894  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
3895  return;
3896 }
3897 
3898 namespace
3899 {
3900  /*SaveResults1Line do single line of output for the save results and save line intensity commands */
3901  /* the number of emission lines across one line of printout */
3902  void SaveLineResults::save(
3903  /* 4 char + null string */
3904  const LinSv *line)
3905  {
3906 
3907  DEBUG_ENTRY( "SaveLineResults::save()" );
3908 
3909  /* if LineWidth is changed then change format in write too */
3910 
3911  /* save results in array so that they can be printed when done */
3912  m_lines[ipLine] = line;
3913 
3914  /* now increment the counter and then check if we have filled the line,
3915  * and so should print it */
3916  ++ipLine;
3917  /* do print now if we are in column mode (one line per line) or if we have filled up
3918  * the line */
3919  if( ( strcmp(::save.chPunRltType,"column") == 0 ) || ipLine == LINEWIDTH )
3920  {
3921  /* "array " is usual array 6 wide, "column" is one line per line */
3922  flush();
3923  }
3924  }
3925 }
3926 
3927 /*SaveGaunts called by save gaunts command to output Gaunt factors */
3928 STATIC void SaveGaunts(FILE* ioPUN)
3929 {
3930  DEBUG_ENTRY( "SaveGaunts()" );
3931 
3932  static const int NENR_GAUNT = 33;
3933  static const int NTE_GAUNT = 20;
3934 
3935  double ener[NENR_GAUNT], ste[NTE_GAUNT], g[NENR_GAUNT][NTE_GAUNT];
3936 
3937  /* this routine is called from the PUNCH GAUNTS command
3938  * it drives the Gaunt factor routine to save gaunts over full range */
3939 
3940  for( int i=0; i < NTE_GAUNT; i++ )
3941  ste[i] = 0.5f*(i+1);
3942 
3943  for( int i=0; i < NENR_GAUNT; i++ )
3944  ener[i] = 0.5f*i - 9.f;
3945 
3946  for( int charge=1; charge <= LIMELM; charge++ )
3947  {
3948  /* energy is log of energy */
3949  for( int ite=0; ite < NTE_GAUNT; ite++ )
3950  {
3951  for( int j=0; j < NENR_GAUNT; j++ )
3952  {
3953  g[j][ite] = t_gaunt::Inst().gauntff( charge, exp10(ste[ite]), exp10(ener[j]) );
3954  }
3955  }
3956 
3957  /* now save out the results */
3958  fprintf( ioPUN, "\tlg(nu)\\lg(Te)" );
3959  for( int i=1; i <= NTE_GAUNT; i++ )
3960  fprintf( ioPUN, "\t%.3e", ste[i-1] );
3961  fprintf( ioPUN, "\n" );
3962 
3963  for( int j=0; j < NENR_GAUNT; j++ )
3964  {
3965  fprintf( ioPUN, "\t%10.3e", ener[j] );
3966  for( int ite=0; ite < NTE_GAUNT; ite++ )
3967  fprintf( ioPUN, "\t%.3e", g[j][ite] );
3968  fprintf( ioPUN, "\n" );
3969  }
3970 
3971  fprintf( ioPUN, "\tlg(nu)/lg(Te)" );
3972  for( int i=0; i < NTE_GAUNT; i++ )
3973  fprintf( ioPUN, "\t%.3e", ste[i] );
3974  fprintf( ioPUN, "\n\n" );
3975 
3976  fprintf( ioPUN, "Below is log(gamma^2), log(u), gff\n" );
3977  /* print log(gamma2), log(u) instead of temp and energy. */
3978 
3979  double z = log10((double)charge);
3980 
3981  for( int i=0; i < NTE_GAUNT; i++ )
3982  {
3983  for( int j=0; j < NENR_GAUNT; j++ )
3984  {
3985  fprintf( ioPUN, "\t%10.3e\t%10.3e\t%10.3e\n",
3986  2.*z + log10(TE1RYD) - ste[i],
3987  log10(TE1RYD) + ener[j] - ste[i],
3988  g[j][i] );
3989  }
3990  }
3991  fprintf( ioPUN, "end of charge = %i\n", charge );
3992  fprintf( ioPUN, "****************************\n" );
3993  }
3994 }
3995 
3996 void SaveGrid(FILE* pnunit, exit_type status)
3997 {
3998  DEBUG_ENTRY( "SaveGrid()" );
3999 
4000  bool lgCreate = ( pnunit == NULL );
4001 
4002  if( lgCreate )
4003  {
4004  // normally the save file will already be open, but it is possible that
4005  // there was an early exit during parsing before the file was opened,
4006  // in which case we open the file here so we can report the error.
4007  string fnam = GridPointPrefix( optimize.nOptimiz ) + save.chFileName[save.ipSaveGrid];
4008  pnunit = open_data( fnam.c_str(), "w" );
4009  }
4010 
4011  if( optimize.nOptimiz == 0 )
4012  {
4013  /* start of line gives abort and warning summary */
4014  fprintf( pnunit, "#Index\tFailure?\tWarnings?\tExit code\t#rank\t#seq" );
4015  /* print start of each variable command line */
4016  for( int i=0; i < grid.nintparm; i++ )
4017  {
4018  string chStr( optimize.chVarFmt[i] );
4019  fprintf( pnunit, "\t%s", chStr.substr(0,9).c_str() );
4020  }
4021  fprintf( pnunit, "\tgrid parameter string\n" );
4022  }
4023  /* abort / warning summary for this sim */
4024  bool lgNoFailure = ( status == ES_SUCCESS || status == ES_WARNINGS );
4025  fprintf( pnunit, "%9.9ld\t%c\t%c\t%20s\t%ld\t%ld",
4026  optimize.nOptimiz,
4027  TorF(!lgNoFailure),
4028  TorF(warnings.lgWarngs),
4029  cpu.i().chExitStatus(status).c_str(),
4030  cpu.i().nRANK(),
4031  grid.seqNum );
4032  /* the grid parameters */
4033  ostringstream chGridParam;
4034  for( int j=0; j < grid.nintparm; j++ )
4035  {
4036  if( j > 0 )
4037  chGridParam << ", ";
4038  chGridParam << fixed << grid.interpParameters[optimize.nOptimiz][j];
4039  fprintf( pnunit, "\t%f", grid.interpParameters[optimize.nOptimiz][j] );
4040  }
4041  fprintf( pnunit, "\t%s\n", chGridParam.str().c_str() );
4042 
4043  if( lgCreate )
4044  fclose( pnunit );
4045 }
SaveLineStuff
STATIC void SaveLineStuff(FILE *ioPUN, const char *chJob, realnum xLimit)
Definition: save_do.cpp:3596
t_cpu_i::nRANK
long nRANK() const
Definition: cpu.h:395
AGN_He1_CS
void AGN_He1_CS(FILE *ioPun)
Definition: iso_solve.cpp:503
monitor_results.h
t_save::lgPunLstIter
bool lgPunLstIter[LIMPUN]
Definition: save.h:389
t_conv::getCounterZone
long getCounterZone(const long type) const
Definition: conv.h:338
t_secondaries::x12tot
realnum x12tot
Definition: secondaries.h:65
t_pressure::PresTotlInit
double PresTotlInit
Definition: pressure.h:52
t_mesh::emm
double emm() const
Definition: mesh.h:93
t_rfield::ConSourceFcnLocal
realnum ** ConSourceFcnLocal
Definition: rfield.h:142
HydroRecCool
double HydroRecCool(long int n, long int ipZ)
Definition: hydroreccool.cpp:10
mole_global
t_mole_global mole_global
Definition: mole.cpp:7
t_rfield::fine_opt_depth
realnum * fine_opt_depth
Definition: rfield.h:393
t_rfield::fine_anu
realnum * fine_anu
Definition: rfield.h:395
t_conv::hist_pres_density
vector< double > hist_pres_density
Definition: conv.h:291
TexcLine
double TexcLine(const TransitionProxy &t)
Definition: transition.cpp:207
CT_OUTW_DIFF
Definition: save_do.cpp:70
t_hmi::H2_Solomon_dissoc_rate_used_H2g
double H2_Solomon_dissoc_rate_used_H2g
Definition: hmi.h:103
t_radius::Radius
double Radius
Definition: radius.h:31
t_radius::depth
double depth
Definition: radius.h:31
t_save::punarg
realnum punarg[LIMPUN][3]
Definition: save.h:372
GrainVar::dstab
vector< double > dstab
Definition: grainvar.h:525
t_thermal::htot
double htot
Definition: thermal.h:169
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:44
t_rfield::line_count
long int * line_count
Definition: rfield.h:59
open_data
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition: cpu.cpp:765
t_opac::ipElement
long int ipElement[LIMELM][LIMELM][7][3]
Definition: opacity.h:223
Wind::dvdr
realnum dvdr
Definition: wind.h:21
PrettyTransmission
double PrettyTransmission(long j, double transmission)
Definition: save_do.cpp:81
Wind::AccelLine
realnum AccelLine
Definition: wind.h:61
thermal
t_thermal thermal
Definition: thermal.cpp:6
prt_LineLabels
void prt_LineLabels(FILE *ioOUT, bool lgPrintAll)
Definition: prt.cpp:112
SaveSpecial
void SaveSpecial(FILE *io, const char *chTime)
Definition: save_special.cpp:12
t_LineSave::SortWL
vector< size_t > SortWL
Definition: lines.h:134
SaveHeat
void SaveHeat(FILE *io)
Definition: heat_save.cpp:18
chIonLbl
string chIonLbl(const TransitionProxy &t)
Definition: transition.cpp:233
ipMAGNESIUM
const int ipMAGNESIUM
Definition: cddefines.h:360
CoolSave
void CoolSave(FILE *io, const char chJob[])
Definition: cool_save.cpp:16
TransitionList::size
size_t size(void) const
Definition: transition.h:331
colden
t_colden colden
Definition: colden.cpp:5
is_odd
bool is_odd(int j)
Definition: cddefines.h:753
t_opac::opacity_abs
double * opacity_abs
Definition: opacity.h:104
t_pressure::PresInteg
realnum PresInteg
Definition: pressure.h:69
ipPun
STATIC long int ipPun
Definition: save_do.cpp:367
null_mole
molecule * null_mole
Definition: mole_species.cpp:50
t_iso_sp::st
qList st
Definition: iso.h:482
t_dynamics::Cool
double Cool()
Definition: dynamics.cpp:2207
t_phycon::te03
double te03
Definition: phycon.h:58
UTALines
TransitionList UTALines("UTALines",&AnonStates)
t_opac::albedo
double * albedo
Definition: opacity.h:113
exp10
double exp10(double x)
Definition: cddefines.h:1368
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:65
t_conv::hist_pres_error
vector< double > hist_pres_error
Definition: conv.h:291
chLineLbl
string chLineLbl(const TransitionProxy &t)
Definition: transition.h:599
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
t_mesh::widflx
double widflx(size_t i) const
Definition: mesh.h:156
t_dense::EdenHCorr
double EdenHCorr
Definition: dense.h:227
input
t_input input
Definition: input.cpp:12
t_radius::PI4_rinner_sq
double PI4_rinner_sq
Definition: radius.h:31
t_grid::lgGrid
bool lgGrid
Definition: grid.h:41
opac
t_opac opac
Definition: opacity.cpp:5
t_save::ResolutionAbs
realnum ResolutionAbs
Definition: save.h:496
t_mole_global::num_calc
int num_calc
Definition: mole.h:362
t_rfield::flux
realnum ** flux
Definition: rfield.h:68
struc
t_struc struc
Definition: struc.cpp:6
t_rfield::DiffuseLineEmission
realnum * DiffuseLineEmission
Definition: rfield.h:193
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
t_iso_sp::RadRec_caseB
double RadRec_caseB
Definition: iso.h:544
hyperfine
t_hyperfine hyperfine
Definition: hyperfine.cpp:5
t_ADfA::ph1
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat_adfa.h:61
t_rfield::chContLabel
vector< string > chContLabel
Definition: rfield.h:213
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:246
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
t_cpu::i
t_cpu_i & i()
Definition: cpu.h:419
NISO
const int NISO
Definition: cddefines.h:311
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:49
t_hcmap::MapZone
long int MapZone
Definition: hcmap.h:20
t_mesh::mesh_md5sum
string mesh_md5sum() const
Definition: mesh.h:112
t_phycon::EnergyIonization
double EnergyIonization
Definition: phycon.h:41
GrainVar::GasCoolColl
double GasCoolColl
Definition: grainvar.h:544
t_rfield::outlin_noplot
realnum * outlin_noplot
Definition: rfield.h:189
save_line
void save_line(FILE *ip, const char *chDo, bool lgEmergent, long ipPun)
Definition: save_line.cpp:116
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:48
t_pressure::PresIntegElecThin
realnum PresIntegElecThin
Definition: pressure.h:75
t_pressure::PresRamCurr
double PresRamCurr
Definition: pressure.h:39
VERSION_TRNCON
static const long VERSION_TRNCON
Definition: save.h:16
cdCautions
void cdCautions(FILE *ioOUT)
Definition: cddrive.cpp:219
t_rfield::flux_total_incident
realnum ** flux_total_incident
Definition: rfield.h:199
TorF
char TorF(bool l)
Definition: cddefines.h:749
DynaSave
void DynaSave(FILE *ipPnunit, char chJob)
Definition: dynamics.cpp:2168
t_pressure::IntegRhoGravity
double IntegRhoGravity
Definition: pressure.h:83
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:46
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:356
t_conv::lgConvPres
bool lgConvPres
Definition: conv.h:192
magnetic
t_magnetic magnetic
Definition: magnetic.cpp:17
t_rfield::ipMaxBolt
long int ipMaxBolt
Definition: rfield.h:230
t_dynamics::lgTimeDependentStatic
bool lgTimeDependentStatic
Definition: dynamics.h:102
EmissionProxy::TauTot
realnum & TauTot() const
Definition: emission.h:490
warnings
t_warnings warnings
Definition: warnings.cpp:11
cdWarnings
void cdWarnings(FILE *ioPNT)
Definition: cddrive.cpp:192
cont_type
cont_type
Definition: save_do.cpp:67
CHIANTI_Upsilon
double CHIANTI_Upsilon(long, long, long, long, double)
Definition: species2.cpp:884
conv
t_conv conv
Definition: conv.cpp:5
SaveGaunts
STATIC void SaveGaunts(FILE *ioPUN)
Definition: save_do.cpp:3928
CT_REFL_DIFF
Definition: save_do.cpp:71
t_optimize::nOptimiz
long int nOptimiz
Definition: optimize.h:250
HFLines
TransitionList HFLines("HFLines",&AnonStates)
t_mole_local::findrk
double findrk(const char buf[]) const
Definition: mole_reactions.cpp:4148
phycon
t_phycon phycon
Definition: phycon.cpp:6
t_phycon::EnthalpyDensity
double EnthalpyDensity
Definition: phycon.h:50
LineSave
t_LineSave LineSave
Definition: lines.cpp:10
t_DoppVel::TurbVel
realnum TurbVel
Definition: doppvel.h:21
SaveSpeciesBands
void SaveSpeciesBands(const long ipPun, const string &speciesLabel, const string &fileBands)
Definition: species_pseudo_cont.cpp:930
t_save::chHashString
char chHashString[INPUT_LINE_LENGTH]
Definition: save.h:416
ipHe2p1P
const int ipHe2p1P
Definition: iso.h:51
SaveLineData
NORETURN void SaveLineData(FILE *io)
Definition: save_linedata.cpp:20
matchGeneric
vector< genericState > matchGeneric(const char *chLabel, bool lgValidate)
Definition: generic_state.cpp:224
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:999
ipRecNetEsc
const int ipRecNetEsc
Definition: cddefines.h:331
ion_recombAGN
void ion_recombAGN(FILE *io)
Definition: ion_recomb.cpp:216
cdSPEC2
void cdSPEC2(int Option, realnum ReturnedSpectrum[])
Definition: save_do.cpp:148
ipFineCont
long ipFineCont(double energy_ryd)
Definition: cont_ipoint.cpp:161
t_cpu_i::big_endian
bool big_endian() const
Definition: cpu.h:354
t_optimize::chVarFmt
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:267
t_geometry::covgeo
realnum covgeo
Definition: geometry.h:45
t_save::lgFLUSH
bool lgFLUSH
Definition: save.h:419
t_magnetic::pressure
double pressure
Definition: magnetic.h:41
t_rfield::outlin
realnum ** outlin
Definition: rfield.h:189
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
t_save::chPunRltType
char chPunRltType[7]
Definition: save.h:444
t_opac::opacity_sct
double * opacity_sct
Definition: opacity.h:107
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:831
PrtMeanIon
void PrtMeanIon(char chType, bool lgDensity, FILE *)
Definition: prt_meanion.cpp:11
t_geometry::FillFac
realnum FillFac
Definition: geometry.h:29
nzone
long int nzone
Definition: cddefines.cpp:14
CT_REFL_INCI
Definition: save_do.cpp:69
t_hmi::HeatH2Dexc_used
double HeatH2Dexc_used
Definition: hmi.h:140
SaveParams::ipPnunit
FILE * ipPnunit
Definition: save.h:195
t_opac::TauScatFace
realnum * TauScatFace
Definition: opacity.h:100
t_called::lgTalk
bool lgTalk
Definition: called.h:12
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:43
DoppVel
t_DoppVel DoppVel
Definition: doppvel.cpp:5
dynamics
t_dynamics dynamics
Definition: dynamics.cpp:42
t_save::lgRealSave
bool lgRealSave[LIMPUN]
Definition: save.h:303
t_phycon::TeProp
double TeProp
Definition: phycon.h:99
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:481
TauLine2
TransitionList TauLine2("TauLine2",&AnonStates)
molecule
Definition: mole.h:142
exit_type
exit_type
Definition: cddefines.h:142
PrtLinePres
void PrtLinePres(FILE *ioPRESSURE)
Definition: prt_linepres.cpp:16
MIN2
#define MIN2(a, b)
Definition: cddefines.h:803
t_pressure::PresTotlCurr
double PresTotlCurr
Definition: pressure.h:46
t_mesh::anu
double anu(size_t i) const
Definition: mesh.h:120
nSpecies
long int nSpecies
Definition: taulines.cpp:22
TempChange
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:32
t_save::lgEmergent
bool lgEmergent[LIMPUN]
Definition: save.h:306
t_LineSave::lines
vector< LinSv > lines
Definition: lines.h:132
ipSULPHUR
const int ipSULPHUR
Definition: cddefines.h:364
t_dynamics::lgStatic_completed
bool lgStatic_completed
Definition: dynamics.h:111
t_save::nsave
long int nsave
Definition: save.h:318
mole_save
void mole_save(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, bool lgCoef, double depth)
Definition: save_species.cpp:527
t_pressure::pres_radiation_lines_curr
double pres_radiation_lines_curr
Definition: pressure.h:61
dense
t_dense dense
Definition: global.cpp:15
GrainVar::GraphiteEmission
vector< realnum > GraphiteEmission
Definition: grainvar.h:579
Singleton< t_yield >::Inst
static t_yield & Inst()
Definition: cddefines.h:209
t_rfield::comup
double * comup
Definition: rfield.h:236
SaveDo
void SaveDo(const char *chTime)
Definition: save_do.cpp:460
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
t_rfield::nflux_with_check
long int nflux_with_check
Definition: rfield.h:49
GrainVar::dstsc
vector< double > dstsc
Definition: grainvar.h:526
atmdat_gaunt.h
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:391
t_dynamics::Heat
double Heat()
Definition: dynamics.cpp:2193
lgSaveOpticalDepths
static bool lgSaveOpticalDepths
Definition: save_do.cpp:3593
RT_diffuse
void RT_diffuse(void)
Definition: rt_diffuse.cpp:35
t_save::lgCumulative
bool lgCumulative[LIMPUN]
Definition: save.h:309
t_rfield::getCoarseTransCoef
const realnum * getCoarseTransCoef()
Definition: rfield.cpp:56
t_radius::PI4_Radius_sq
double PI4_Radius_sq
Definition: radius.h:31
t_ionbal::RateRecomTot
double ** RateRecomTot
Definition: ionbal.h:194
t_grid::interpParameters
realnum ** interpParameters
Definition: grid.h:31
t_magnetic::EnthalpyDensity
double EnthalpyDensity
Definition: magnetic.h:38
CT_GRN_TOT
Definition: save_do.cpp:76
t_save::optname
string optname[LIMPUN]
Definition: save.h:375
wind
Wind wind
Definition: wind.cpp:5
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
t_geometry::lgSphere
bool lgSphere
Definition: geometry.h:34
t_conv::hist_temp_cool
vector< double > hist_temp_cool
Definition: conv.h:296
t_grid::nintparm
long nintparm
Definition: grid.h:58
iteration
long int iteration
Definition: cddefines.cpp:16
t_save::SpeciesBandFile
string SpeciesBandFile[LIMPUN]
Definition: save.h:513
t_rfield::otslin
realnum * otslin
Definition: rfield.h:183
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:854
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:554
t_save::LinEvery
long int LinEvery
Definition: save.h:480
SaveSpeciesOptDep
void SaveSpeciesOptDep(const long int ipPun, const string &speciesLabel)
Definition: save_species.cpp:441
t_radius::drad
double drad
Definition: radius.h:31
save_opacity
void save_opacity(FILE *io, long int np)
Definition: save_opacity.cpp:23
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:8
abund
t_abund abund
Definition: abund.cpp:5
GrainVar::TotalEden
double TotalEden
Definition: grainvar.h:529
Wind::AccelGravity
realnum AccelGravity
Definition: wind.h:49
t_save::lg_separate_iterations
bool lg_separate_iterations[LIMPUN]
Definition: save.h:334
t_rfield::ConEmitLocal
realnum ** ConEmitLocal
Definition: rfield.h:139
t_pressure::pinzon
realnum pinzon
Definition: pressure.h:69
geometry
t_geometry geometry
Definition: geometry.cpp:5
t_conv::getCounterName
const char * getCounterName(const long type) const
Definition: conv.h:342
Wind::AccelCont
realnum AccelCont
Definition: wind.h:55
ipIRON
const int ipIRON
Definition: cddefines.h:374
colliders
ColliderList colliders(dense)
ipoint
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:15
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:570
t_hcmap::lgMapDone
bool lgMapDone
Definition: hcmap.h:36
t_dense::pden
realnum pden
Definition: dense.h:108
t_yield::nelec_eject
long nelec_eject(long n, long i, long ns) const
Definition: yield.h:55
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:477
SaveGrid
void SaveGrid(FILE *pnunit, exit_type status)
Definition: save_do.cpp:3996
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
GrainVar::GrainEmission
vector< realnum > GrainEmission
Definition: grainvar.h:578
t_yield::elec_eject_frac
realnum elec_eject_frac(long n, long i, long ns, long ne) const
Definition: yield.h:48
POW2
#define POW2
Definition: cddefines.h:979
CT_GRN_SIL
Definition: save_do.cpp:74
t_isoCTRL::nLyman
long int nLyman[NISO]
Definition: iso.h:352
energy
double energy(const genericState &gs)
Definition: generic_state.cpp:49
ipH1s
const int ipH1s
Definition: iso.h:29
t_hmi::HeatH2Dish_used
double HeatH2Dish_used
Definition: hmi.h:140
t_conv::nPres2Ioniz
long int nPres2Ioniz
Definition: conv.h:145
STATIC
#define STATIC
Definition: cddefines.h:118
t_save::chSaveArgs
char chSaveArgs[LIMPUN][5]
Definition: save.h:383
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:610
flxCell
double flxCell(long j, long nEmType, cont_type ct, bool lgForceConserve=false, bool lgPrtIsotropicCont=true, const realnum *trans_coef_total=NULL)
Definition: save_do.cpp:94
diatomics::lgEnabled
bool lgEnabled
Definition: h2_priv.h:352
EmissionProxy::Ploss
realnum Ploss() const
Definition: emission.h:130
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:447
ipExtraLymanLines
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
Save_Line_RT
void Save_Line_RT(FILE *ip)
Definition: save_line.cpp:334
t_save::ncSaveSkip
long int ncSaveSkip
Definition: save.h:484
t_colden::H0_21cm_lower
double H0_21cm_lower
Definition: colden.h:67
FindStrongestLineLabels
STATIC void FindStrongestLineLabels(void)
Definition: save_do.cpp:273
mole
t_mole_local mole
Definition: mole.cpp:8
t_conv::hist_pres_current
vector< double > hist_pres_current
Definition: conv.h:291
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:789
t_dense::EdenTrue
double EdenTrue
Definition: dense.h:232
saveFITSfile
void saveFITSfile(FILE *io, int option, realnum Elo=0.f, realnum Ehi=0.f, realnum Enorm=0.f)
Definition: save_fits.cpp:86
pressure
t_pressure pressure
Definition: pressure.cpp:9
rfield
t_rfield rfield
Definition: rfield.cpp:9
mean
t_mean mean
Definition: mean.cpp:16
TransitionProxy::ipCont
long & ipCont() const
Definition: transition.h:489
t_conv::hist_temp_temp
vector< double > hist_temp_temp
Definition: conv.h:296
t_rfield::ConInterOut
realnum * ConInterOut
Definition: rfield.h:154
realnum
float realnum
Definition: cddefines.h:124
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:468
NUM_OUTPUT_TYPES
const int NUM_OUTPUT_TYPES
Definition: grid.h:22
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
t_colden::H0_21cm_upper
double H0_21cm_upper
Definition: colden.h:66
Wind::AccelTotalOutward
realnum AccelTotalOutward
Definition: wind.h:52
t_input::echo
void echo(FILE *ipOUT)
Definition: input.cpp:199
t_colden::H0_ov_Tspin
double H0_ov_Tspin
Definition: colden.h:52
EmissionProxy::Pesc
realnum & Pesc() const
Definition: emission.h:580
t_save::nSaveEveryZone
long int nSaveEveryZone[LIMPUN]
Definition: save.h:380
CT_OUTW_LIN
Definition: save_do.cpp:72
max
long max(int a, long b)
Definition: cddefines.h:817
t_save::ipSaveGrid
long ipSaveGrid
Definition: save.h:284
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:435
t_opac::telec
realnum telec
Definition: opacity.h:188
t_rfield::otscon
realnum * otscon
Definition: rfield.h:183
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:482
t_save::lgLineListRatio
bool lgLineListRatio[LIMPUN]
Definition: save.h:265
molecule::index
int index
Definition: mole.h:194
t_cpu_i::lgMPI_talk
bool lgMPI_talk() const
Definition: cpu.h:397
t_timesc::time_H2_Dest_here
double time_H2_Dest_here
Definition: timesc.h:48
ipRecRad
const int ipRecRad
Definition: cddefines.h:333
t_save::chLineListLabel
vector< string > chLineListLabel[LIMPUN]
Definition: save.h:261
t_radius::lgRadiusKnown
bool lgRadiusKnown
Definition: radius.h:122
t_hcmap::lgMapBeingDone
bool lgMapBeingDone
Definition: hcmap.h:33
t_opac::TauAbsFace
realnum * TauAbsFace
Definition: opacity.h:100
t_mesh::anu2
double anu2(size_t i) const
Definition: mesh.h:124
h2
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
CT_OUTW_INCI
Definition: save_do.cpp:68
GetDopplerWidth
realnum GetDopplerWidth(realnum massAMU)
Definition: temp_change.cpp:376
safe_div
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1011
t_rfield::OccNumbIncidCont
realnum * OccNumbIncidCont
Definition: rfield.h:117
t_radius::depth_mid_zone
double depth_mid_zone
Definition: radius.h:31
iterations
t_iterations iterations
Definition: iterations.cpp:6
ipRecEsc
const int ipRecEsc
Definition: cddefines.h:329
column
double column(const genericState &gs)
Definition: generic_state.cpp:15
map_do
void map_do(FILE *io, const char *chType)
Definition: hcmap.cpp:25
t_pressure::PresGasCurr
double PresGasCurr
Definition: pressure.h:46
nWindLine
long nWindLine
Definition: cdinit.cpp:19
t_pressure::PresTotlError
double PresTotlError
Definition: pressure.h:46
optimize
t_optimize optimize
Definition: optimize.cpp:6
PrtLogLin
double PrtLogLin(double value)
Definition: save_do.cpp:370
lgMustPrintHeader
static bool lgMustPrintHeader
Definition: save_line.cpp:270
t_rfield::chLineLabel
vector< string > chLineLabel
Definition: rfield.h:210
grid
t_grid grid
Definition: grid.cpp:5
cdExecTime
double cdExecTime()
Definition: cddrive.cpp:478
radius
t_radius radius
Definition: radius.cpp:5
t_thermal::heating
double heating(long nelem, long ion)
Definition: thermal.h:186
timesc
t_timesc timesc
Definition: timesc.cpp:7
t_hmi::H2_photodissoc_used_H2g
double H2_photodissoc_used_H2g
Definition: hmi.h:119
AnuUnit
double AnuUnit(realnum energy)
Definition: service.cpp:197
t_mesh::anumin
double anumin(size_t i) const
Definition: mesh.h:148
t_save::params
SaveParams params[LIMPUN]
Definition: save.h:278
t_rfield::reflin
realnum ** reflin
Definition: rfield.h:196
t_ionbal::PhotoRate_Shell
double **** PhotoRate_Shell
Definition: ionbal.h:112
t_rfield::ConEmitOut
realnum ** ConEmitOut
Definition: rfield.h:151
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:160
t_gaunt::gauntff
double gauntff(long Z, double Te, double anu)
Definition: atmdat_gaunt.cpp:135
t_rfield::fine_opac_zone
realnum * fine_opac_zone
Definition: rfield.h:391
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_rfield::extin_mag_V_point
double extin_mag_V_point
Definition: rfield.h:258
chDummy
char * chDummy
Definition: save_do.cpp:458
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:76
t_hmi::UV_Cont_rel2_Draine_DB96_depth
realnum UV_Cont_rel2_Draine_DB96_depth
Definition: hmi.h:84
t_radius::Conv2PrtInten
double Conv2PrtInten
Definition: radius.h:153
ipSILICON
const int ipSILICON
Definition: cddefines.h:362
t_radius::Radius_mid_zone
double Radius_mid_zone
Definition: radius.h:31
t_rfield::ipV_filter
long int ipV_filter
Definition: rfield.h:240
ipH2p
const int ipH2p
Definition: iso.h:31
t_save::lgSaveEveryZone
bool lgSaveEveryZone[LIMPUN]
Definition: save.h:379
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
EmissionProxy::Pdest
realnum & Pdest() const
Definition: emission.h:600
t_save::lgLinEvery
bool lgLinEvery
Definition: save.h:481
t_save::SaveHeaderDone
void SaveHeaderDone(int ipPun)
Definition: save.h:364
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:613
sprt_wl
void sprt_wl(char *chString, realnum wl)
Definition: prt.cpp:56
t_save::lgPrtIsotropicCont
bool lgPrtIsotropicCont[LIMPUN]
Definition: save.h:315
t_timesc::sound_speed_adiabatic
double sound_speed_adiabatic
Definition: timesc.h:58
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:431
t_iterations::lgLastIt
bool lgLastIt
Definition: iterations.h:47
t_save::lgFITS
bool lgFITS[LIMPUN]
Definition: save.h:392
OccupationNumberLine
double OccupationNumberLine(const TransitionProxy &t)
Definition: transition.cpp:177
t_hyperfine::Tspin21cm
double Tspin21cm
Definition: hyperfine.h:56
EmissionProxy::TauCon
realnum & TauCon() const
Definition: emission.h:510
ipH2s
const int ipH2s
Definition: iso.h:30
t_conv::lgConvTemp
bool lgConvTemp
Definition: conv.h:189
ipALUMINIUM
const int ipALUMINIUM
Definition: cddefines.h:361
mole_dominant_rates
void mole_dominant_rates(const vector< const molecule * > &debug_list, FILE *ioOut, bool lgPrintReagents, size_t NPRINT, double fprint)
Definition: save_species.cpp:612
t_phycon::EnergyBinding
double EnergyBinding
Definition: phycon.h:54
t_dynamics::Rate
double Rate
Definition: dynamics.h:77
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1067
t_save::wlLineList
vector< realnum > wlLineList[LIMPUN]
Definition: save.h:263
ChargTranPun
void ChargTranPun(FILE *ipPnunit, char *chSave)
Definition: atmdat_char_tran.cpp:1728
density
double density(const genericState &gs)
Definition: generic_state.cpp:27
t_rfield::extin_mag_V_extended
double extin_mag_V_extended
Definition: rfield.h:262
t_save::lgHashEndIter
bool lgHashEndIter[LIMPUN]
Definition: save.h:412
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:64
ExtraLymanLines
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
LIMELM
const int LIMELM
Definition: cddefines.h:308
cdColm
int cdColm(const char *chLabel, long int ion, double *theocl)
Definition: cddrive.cpp:592
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:50
t_rfield::nfine
long nfine
Definition: rfield.h:387
t_radius::drad_x_fillfac
double drad_x_fillfac
Definition: radius.h:77
Energy::Ryd
double Ryd() const
Definition: energy.h:33
molezone::den
double den
Definition: mole.h:421
t_save::lgPrtOldStyleLogs
bool lgPrtOldStyleLogs[LIMPUN]
Definition: save.h:290
TransitionProxy::Coll
CollisionProxy Coll() const
Definition: transition.h:463
diatomics::H2_PunchLineStuff
void H2_PunchLineStuff(FILE *io, realnum xLimit, long index)
Definition: mole_h2_io.cpp:1158
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:350
generic_state.h
t_rfield::comdn
double * comdn
Definition: rfield.h:236
t_phycon::te10
double te10
Definition: phycon.h:58
t_dense::xMassDensity
realnum xMassDensity
Definition: dense.h:101
diatomics::H2_PunchDo
void H2_PunchDo(FILE *io, char chJOB[], const char chTime[], long int ipPun)
Definition: mole_h2_io.cpp:1195
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:583
t_hmi::H2_total
double H2_total
Definition: hmi.h:25
t_save::whichDiatomToPrint
diatomics * whichDiatomToPrint[LIMPUN]
Definition: save.h:322
ConvPresTempEdenIoniz
int ConvPresTempEdenIoniz(void)
Definition: conv_pres_temp_eden_ioniz.cpp:23
t_mesh::egamry
double egamry() const
Definition: mesh.h:97
t_conv::hist_temp_heat
vector< double > hist_temp_heat
Definition: conv.h:296
t_phycon::TeInit
double TeInit
Definition: phycon.h:99
t_save::ipEmisFreq
long ipEmisFreq[LIMPUN]
Definition: save.h:502
t_grid::lgOutputTypeOn
bool lgOutputTypeOn[NUM_OUTPUT_TYPES]
Definition: grid.h:66
dBaseSpecies
vector< species > dBaseSpecies
Definition: taulines.cpp:15
t_dense::eden
double eden
Definition: dense.h:201
t_opac::TauAbsGeo
realnum ** TauAbsGeo
Definition: opacity.h:91
SaveSpecies
void SaveSpecies(FILE *ioPUN, long int ipPun)
Definition: save_species.cpp:35
TransitionProxy::EnergyRyd
double EnergyRyd() const
Definition: transition.h:95
t_phycon::EnergyExcitation
double EnergyExcitation
Definition: phycon.h:47
NCHLAB
const int NCHLAB
Definition: cddefines.h:304
MAX2
#define MAX2(a, b)
Definition: cddefines.h:824
t_ionbal::RateIonizTot
double RateIonizTot(long nelem, long ion) const
Definition: ionbal.cpp:224
t_hmi::UV_Cont_rel2_Habing_TH85_depth
realnum UV_Cont_rel2_Habing_TH85_depth
Definition: hmi.h:74
t_opac::ExpmTau
realnum * ExpmTau
Definition: opacity.h:145
lgCheckMonitors
bool lgCheckMonitors(FILE *ioMONITOR)
Definition: monitor_results.cpp:1703
EmissionProxy::damp
realnum & damp() const
Definition: emission.h:620
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
t_phycon::te70
double te70
Definition: phycon.h:58
t_mean::xIonMean
multi_arr< double, 4 > xIonMean
Definition: mean.h:17
t_Heavy::chShell
char chShell[7][3]
Definition: heavy.h:31
t_save::chSpeciesDominantRates
string chSpeciesDominantRates[LIMPUN]
Definition: save.h:498
t_save::chConSavEnr
const char * chConSavEnr[LIMPUN]
Definition: save.h:405
t_mole_global::list
MoleculeList list
Definition: mole.h:365
t_DoppVel::lgTurb_pressure
bool lgTurb_pressure
Definition: doppvel.h:42
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:33
prt_line_inlist
void prt_line_inlist(FILE *ioOUT, const char *label, realnum wvlng)
Definition: prt.cpp:173
t_dynamics::Source
double ** Source
Definition: dynamics.h:80
AGN_Hemis
STATIC void AGN_Hemis(FILE *ioPUN)
Definition: save_do.cpp:3769
SaveSpeciesPseudoCont
void SaveSpeciesPseudoCont(const long ipPun, const string &speciesLabel)
Definition: species_pseudo_cont.cpp:406
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
t_save::chSaveSpecies
vector< string > chSaveSpecies[LIMPUN]
Definition: save.h:385
t_dynamics::time_elapsed
double time_elapsed
Definition: dynamics.h:105
t_phycon::telogn
double telogn[7]
Definition: phycon.h:86
t_conv::ntypes
size_t ntypes(void) const
Definition: conv.h:305
ipCARBON
const int ipCARBON
Definition: cddefines.h:354
t_rfield::TotDiff2Pht
realnum * TotDiff2Pht
Definition: rfield.h:177
t_grid::seqNum
long seqNum
Definition: grid.h:69
CT_GRN_GRA
Definition: save_do.cpp:75
DynaPunchTimeDep
void DynaPunchTimeDep(FILE *ipPnunit, const char *chJob)
Definition: dynamics.cpp:2059
save_average
void save_average(long int ipPun)
Definition: save_average.cpp:128
hcmap
t_hcmap hcmap
Definition: hcmap.cpp:23
GammaPrt
void GammaPrt(long int ipLoEnr, long int ipHiEnr, long int ipOpac, FILE *ioFILE, double total, double threshold)
Definition: cont_gammas.cpp:244
SaveLineIntensity
STATIC void SaveLineIntensity(FILE *ioPUN, long int ipPun, realnum Threshold)
Definition: save_do.cpp:3545
t_LineSave::nsum
long int nsum
Definition: lines.h:87
GridPointPrefix
string GridPointPrefix(int n)
Definition: mpi_utilities.h:154
t_pressure::pinzon_PresIntegElecThin
realnum pinzon_PresIntegElecThin
Definition: pressure.h:75
t_struc::testr
realnum * testr
Definition: struc.h:25
CT_REFL_LIN
Definition: save_do.cpp:73
GrainVar::SilicateEmission
vector< realnum > SilicateEmission
Definition: grainvar.h:580
t_mesh::anumax
double anumax(size_t i) const
Definition: mesh.h:152
t_pressure::PresTurbCurr
double PresTurbCurr
Definition: pressure.h:42
fixit
#define fixit(a)
Definition: cddefines.h:417
t_rfield::ConEmitReflec
realnum ** ConEmitReflec
Definition: rfield.h:145
gv
GrainVar gv
Definition: grainvar.cpp:5
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
hmi
t_hmi hmi
Definition: hmi.cpp:5
cpu
static t_cpu cpu
Definition: cpu.h:427
ShowMe
void ShowMe(void)
Definition: service.cpp:205
dBaseTrans
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
save
t_save save
Definition: save.cpp:5
t_phycon::te
double te
Definition: phycon.h:21
t_rfield::chCumuType
char chCumuType[5]
Definition: rfield.h:335
t_thermal::dCooldT
double dCooldT
Definition: thermal.h:139
plankf
double plankf(long int ip)
Definition: service.cpp:1760
t_save::Resolution
realnum Resolution
Definition: save.h:493
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:349
SaveResults
STATIC void SaveResults(FILE *ioPUN)
Definition: save_do.cpp:3838
t_rfield::nflux
long int nflux
Definition: rfield.h:46
t_save::emisfreq
Energy emisfreq[LIMPUN]
Definition: save.h:501
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:690
flux_correct_isotropic
double flux_correct_isotropic(const bool lgSaveIsotr, const int nEmType, const int iflux)
Definition: rfield.cpp:112
t_colden::colden
realnum colden[NCOLD]
Definition: colden.h:32
CT_INCI
Definition: save_do.cpp:67
t_save::lgSaveHeader
bool lgSaveHeader(int ipPun) const
Definition: save.h:360
t_rfield::ConRefIncid
realnum ** ConRefIncid
Definition: rfield.h:157
t_save::ipConPun
long int ipConPun
Definition: save.h:408
LinSv
Definition: lines.h:157
t_save::nLineList
long nLineList[LIMPUN]
Definition: save.h:259
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:529
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:519
t_save::chSave
char chSave[LIMPUN][5]
Definition: save.h:321
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:468
t_mesh::getResolutionScaleFactor
double getResolutionScaleFactor() const
Definition: mesh.h:101
Wind::windv
realnum windv
Definition: wind.h:18
ipLine
static long int * ipLine
Definition: prt_linesum.cpp:14
CollisionProxy::ColUL
double ColUL(const ColliderList &colls) const
Definition: collision.h:106
EmissionProxy::TauIn
realnum & TauIn() const
Definition: emission.h:470
t_warnings::lgWarngs
bool lgWarngs
Definition: warnings.h:44
t_LineSave::ipass
long int ipass
Definition: lines.h:96
called
t_called called
Definition: called.cpp:4
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:363
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
GrainVar::rate_h2_form_grains_used_total
double rate_h2_form_grains_used_total
Definition: grainvar.h:574
t_cpu_i::chExitStatus
const string & chExitStatus(exit_type s) const
Definition: cpu.h:403
ipSODIUM
const int ipSODIUM
Definition: cddefines.h:359
lgPopsFirstCall
static bool lgPopsFirstCall
Definition: save_do.cpp:3593
t_thermal::ctot
double ctot
Definition: thermal.h:130
Wind::fmul
realnum fmul
Definition: wind.h:65
t_save::chFileName
vector< string > chFileName
Definition: save.h:281
Save1Line
void Save1Line(const TransitionProxy &t, FILE *io, realnum xLimit, long index, realnum DopplerWidth)
Definition: save_do.cpp:3711
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:166
EmissionProxy::pump
double & pump() const
Definition: emission.h:530
Generated by   doxygen 1.8.5