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