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
70 STATIC void FindStrongestLineLabels( void )
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  {
715  realnum corr = save.ResolutionAbs*rfield.widflx(j)/rfield.anu(j);
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");
787  cdEXIT(EXIT_FAILURE);
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 )
861  GridGatherInCloudy();
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
871  FindStrongestLineLabels();
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 =
889  ( iterations.lgLastIt || !save.lgPunLstIter[ipPun] ||
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 ) ||
895  ( dynamics.lgTimeDependentStatic && dynamics.lgStatic_completed ) );
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",
906  log10(MAX2(SMALLFLOAT,dense.gas_phase[ipHYDROGEN])) );
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 */
926  radius.depth_mid_zone,
927  hyperfine.Tspin21cm ,
928  phycon.te ,
929  /* temperature from Lya - 21 cm optical depth ratio */
930  3.84e-7* iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon() /
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() ,
935  OccupationNumberLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ),
936  HFLines[0].Emis().TauCon() ,
937  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon(),
938  HFLines[0].Emis().PopOpc(),
939  /* term in () is density (cm-3) of 1s, so this is n(1s) / Ts */
940  (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop())/SDIV( hyperfine.Tspin21cm),
941  /* why was above multiplied by this following term? */
942  /* *HFLines[0].EnergyErg/BOLTZMANN/4.,*/
943  HFLines[0].Emis().TauIn(),
944  TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ) ,
945  colden.H0_ov_Tspin,
946  /*>>chng 27 mar, GS, integrated 21cm spin temperature*/
947  colden.H0_21cm_lower,
948  colden.H0_21cm_upper,
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 */
964  radius.depth_mid_zone,
965  /* cooling timescale */
966  dense.pden*BOLTZMANN*1.5*phycon.te/ thermal.htot,
967  /* H2 destruction timescale */
968  timesc.time_H2_Dest_here,
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();
997  cdEXIT(EXIT_FAILURE);
998  }
999  }
1000  }
1001 
1002  else if( strcmp(save.chSave[ipPun],"MONI") == 0 )
1003  {
1004  if( lgLastOnly )
1005  {
1006  /* save the monitor output */
1007  lgCheckMonitors(save.params[ipPun].ipPnunit);
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",
1151  radius.depth_mid_zone,
1152  phycon.EnthalpyDensity,
1153  phycon.EnergyExcitation,
1154  phycon.EnergyIonization,
1155  phycon.EnergyBinding ,
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) *
1237  rfield.anu2(j)*radius.PI4_rinner_sq*
1238  EN1RYD;
1239 
1240  // a value < 0. indicates that energy should be conserved
1241  realnum resolution = ( save.Resolution < realnum(0.) ) ?
1242  rfield.anu(j)/rfield.widflx(j) : save.Resolution;
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/
1251  rfield.widflx(j)*radius.PI4_Radius_sq *
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*
1280  radius.PI4_Radius_sq*geometry.covgeo);
1281 
1282  /* outward lines */
1283  fprintf( save.params[ipPun].ipPnunit,"%.3e\t", rfield.anu(j)*resolution*rfield.outlin[nEmType][j]*EN1RYD*
1284  radius.PI4_Radius_sq*geometry.covgeo );
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 )
1303  SaveNewContinuum( save.params[ipPun].ipPnunit);
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)*
1335  EN1RYD/rfield.widflx(j) *radius.PI4_Radius_sq*geometry.covgeo;
1336 
1337  double fout = gv.SilicateEmission[j]*rfield.anu2(j)*
1338  EN1RYD/rfield.widflx(j) *radius.PI4_Radius_sq*geometry.covgeo;
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)*
1347  EN1RYD/rfield.widflx(j) *radius.PI4_Radius_sq*geometry.covgeo );
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" );
1364  cdEXIT(EXIT_FAILURE);
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.) ) ?
1371  rfield.anu(j)/rfield.widflx(j) : save.Resolution;
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",
1403  radius.depth_mid_zone,
1404  conv.nPres2Ioniz,
1405  pressure.PresTotlCurr,
1406  pressure.PresTotlError*100.,
1407  dense.EdenTrue,
1408  dense.eden,
1409  (dense.EdenTrue - dense.eden)*100./dense.EdenTrue,
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.) ) ?
1446  rfield.anu(j)/rfield.widflx(j) : save.Resolution;
1447 
1448  double EmisLin = resolution*EN1RYD*
1449  rfield.DiffuseLineEmission[j]*rfield.anu(j);
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.) ) ?
1478  rfield.anu(0)/rfield.widflx(0) : save.Resolution;
1479  double EmisLin = resolution*EN1RYD*
1480  rfield.DiffuseLineEmission[0]*rfield.anu(0);
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.) ) ?
1489  rfield.anu(j)/rfield.widflx(j) : save.Resolution;
1490  double EmisLin = resolution*EN1RYD*
1491  rfield.DiffuseLineEmission[j]*rfield.anu(j);
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.) ) ?
1513  rfield.anu(j)/rfield.widflx(j) : save.Resolution;
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.) ) ?
1734  rfield.anu(j)/rfield.widflx(j) : save.Resolution;
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/
1744  rfield.widflx(j)*radius.PI4_Radius_sq *
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",
1790  radius.Radius_mid_zone,
1791  radius.depth_mid_zone,
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] +
1844  rfield.ConInterOut[j])*opac.opacity_abs[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.) ) ?
1863  rfield.anu(j)/rfield.widflx(j) : save.Resolution;
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,
1868  rfield.ConInterOut[j]/rfield.widflx(j)*rfield.anu2(j)*radius.PI4_Radius_sq*EN1RYD,
1869  resolution*(rfield.outlin[0][j]+rfield.outlin_noplot[j])*rfield.anu(j)*radius.PI4_Radius_sq*EN1RYD,
1870  (rfield.ConInterOut[j]/rfield.widflx(j)*
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] );
1909  u.x = rfield.getResolutionScaleFactor();
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",
1917  (rfield.nflux+save.ncSaveSkip-1)/save.ncSaveSkip );
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)*
1948  rfield.anu2(j)*radius.PI4_Radius_sq*EN1RYD*geometry.covgeo;
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",
1984  radius.depth_mid_zone );
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",
2047  radius.depth_mid_zone );
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",
2079  radius.depth_mid_zone );
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",
2102  radius.depth_mid_zone );
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",
2126  radius.depth_mid_zone );
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",
2151  radius.depth_mid_zone ,
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",
2181  radius.depth_mid_zone );
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",
2205  radius.depth_mid_zone );
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",
2307  radius.Radius, radius.depth ,wind.windv/1e5,
2308  wind.dvdr,
2309  dense.gas_phase[ipHYDROGEN], dense.eden , phycon.te,
2310  wind.AccelLine , wind.AccelCont ,
2311  wind.fmul ,
2312  // acceleration in this zone due to electron scattering,
2313  // if incident SED was not attenuated
2314  pressure.pinzon_PresIntegElecThin/dense.xMassDensity/radius.drad_x_fillfac ,
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] ,
2324  dense.xIonDense[ipHYDROGEN][0] , dense.xIonDense[ipHYDROGEN][1] ,
2325  dense.xIonDense[ipHELIUM][0] , dense.xIonDense[ipHELIUM][1] ,
2326  dense.xIonDense[ipHELIUM][2] ,
2327  dense.xIonDense[ipCARBON][1] , dense.xIonDense[ipCARBON][2] ,
2328  dense.xIonDense[ipCARBON][3] ,
2329  dense.xIonDense[ipOXYGEN][0] , dense.xIonDense[ipOXYGEN][1] ,
2330  dense.xIonDense[ipOXYGEN][2] , dense.xIonDense[ipOXYGEN][3] ,
2331  dense.xIonDense[ipOXYGEN][4] , dense.xIonDense[ipOXYGEN][5] ,
2332  dense.xIonDense[ipOXYGEN][6] , dense.xIonDense[ipOXYGEN][7] ,
2333  mean.xIonMean[0][ipMAGNESIUM][1][0] , dense.xIonDense[ipMAGNESIUM][1]);
2334  }
2335  }
2336 
2337  /* save spectra in fits format */
2338  else if( strcmp(save.chSave[ipPun],"FITS") == 0 )
2339  {
2340  if( lgLastOnly )
2341  {
2342  saveFITSfile( save.params[ipPun].ipPnunit, NUM_OUTPUT_TYPES );
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,
2445  conv.hist_pres_density[k],
2446  conv.hist_pres_current[k],
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 ,
2518  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).WLAng() );
2519  fprintf( save.params[ipPun].ipPnunit, "\t" );
2520  prt_wl( save.params[ipPun].ipPnunit ,
2521  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P1,ipHe1s1S).WLAng() );
2522  fprintf( save.params[ipPun].ipPnunit, "\t" );
2523  prt_wl( save.params[ipPun].ipPnunit ,
2524  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P2,ipHe1s1S).WLAng() );
2525  fprintf( save.params[ipPun].ipPnunit, "\t" );
2526  prt_wl( save.params[ipPun].ipPnunit ,
2527  iso_sp[ipHE_LIKE][nelem].trans(ipHe2s3S,ipHe1s1S).WLAng() );
2528  fprintf( save.params[ipPun].ipPnunit, "\t" );
2529  prt_wl( save.params[ipPun].ipPnunit ,
2530  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P2,ipHe2s3S).WLAng() );
2531  fprintf( save.params[ipPun].ipPnunit, "\t" );
2532  prt_wl( save.params[ipPun].ipPnunit ,
2533  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P1,ipHe2s3S).WLAng() );
2534  fprintf( save.params[ipPun].ipPnunit, "\t" );
2535  prt_wl( save.params[ipPun].ipPnunit ,
2536  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P0,ipHe2s3S).WLAng() );
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()/
2549  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Coll().ColUL( colliders );
2550  fprintf( save.params[ipPun].ipPnunit,
2551  " %.5e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
2552  radius.depth_mid_zone,
2553  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn(),
2554  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop(),
2555  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop(),
2556  phycon.te, iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().damp(), eps );
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",
2569  radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], dense.eden,
2570  dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
2571  dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
2572  hmi.H2_total/dense.gas_phase[ipHYDROGEN],
2573  findspecieslocal("H2+")->den/dense.gas_phase[ipHYDROGEN],
2574  findspecieslocal("H3+")->den/dense.gas_phase[ipHYDROGEN],
2575  findspecieslocal("H-")->den/dense.gas_phase[ipHYDROGEN] );
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 +=
2594  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Aul()*
2595  (iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Ploss());
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,
2614  ionbal.RateRecomTot[ipHYDROGEN][0] );
2615 
2616  fprintf( save.params[ipPun].ipPnunit, "\t%.2e",
2617  iso_sp[ipH_LIKE][ipHYDROGEN].RadRec_caseB );
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",
2644  radius.depth_mid_zone,
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 -
2665  iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local; ++ipHi )
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",
2694  radius.depth_mid_zone,
2695  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn(),
2696  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot(),
2697  popul,
2698  texc,
2699  phycon.te,
2700  texc/phycon.te ,
2701  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pesc(),
2702  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest(),
2703  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().pump(),
2704  opac.opacity_abs[iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont()-1],
2705  opac.albedo[iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont()-1] );
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",
2767  radius.depth_mid_zone,
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 &&
2978  (nzone/save.LinEvery)*save.LinEvery != nzone )
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) ||
2988  (nzone/save.LinEvery)*save.LinEvery == nzone )
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 */
3088  radius.depth_mid_zone,
3089  /* A_V for an extended source */
3090  0.00,
3091  /* A_V for a point source */
3092  rfield.extin_mag_V_point,
3093  /* temperature */
3094  phycon.te ,
3095  dense.xIonDense[ipHYDROGEN][0],
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],
3106  dense.xIonDense[ipHYDROGEN][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,
3118  colden.colden[ipCOL_elec],
3119  col_heii,
3120  findspecieslocal("H+")->column,
3121  findspecieslocal("H3+")->column,
3122  hmi.H2_Solomon_dissoc_rate_used_H2g ,
3123  gv.rate_h2_form_grains_used_total,
3124  hmi.H2_photodissoc_used_H2g,
3125  hmi.UV_Cont_rel2_Draine_DB96_depth,
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 */
3143  hmi.HeatH2Dish_used,
3144  /* H2vH-Heating due to collisions with H2 */
3145  hmi.HeatH2Dexc_used ,
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 */
3151  thermal.heating(ipMAGNESIUM,0),
3152  thermal.heating(ipSULPHUR,0),
3153  thermal.heating(ipSILICON,0),
3154  thermal.heating(ipIRON,0),
3155  thermal.heating(ipSODIUM,0),
3156  thermal.heating(ipALUMINIUM,0),
3157  thermal.heating(ipCARBON,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");
3281  cdEXIT(EXIT_FAILURE);
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)),
3364  opac.TauAbsFace[j]+opac.TauScatFace[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 ),
3469  PrtLogLin(thermal.htot-dynamics.Heat() ) );
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",
3481  PrtLogLin(dense.gas_phase[ipHYDROGEN]),
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",
3490  PrtLogLin(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])),
3491  PrtLogLin(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])) );
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",
3554  radius.depth_mid_zone,
3555  /* total hydrogen column density, all forms */
3556  colden.colden[ipCOL_HTOT],
3557  phycon.te,
3558  /* fraction of H that is atomic */
3559  dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
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 */
3564  dense.xIonDense[ipCARBON][0]/SDIV(dense.gas_phase[ipCARBON]),
3565  findspecieslocal("CO")->den/SDIV(dense.gas_phase[ipCARBON]),
3566  findspecieslocal("H2O")->den/SDIV(dense.gas_phase[ipOXYGEN]),
3567  /* hmi.UV_Cont_rel2_Habing_TH85 is field relative to Habing background, dimensionless */
3568  hmi.UV_Cont_rel2_Habing_TH85_depth);
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",
3627  radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN],
3628  dense.eden, thermal.htot, wind.AccelTotalOutward, geometry.FillFac );
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 */
3640  radius.depth_mid_zone,
3641  /*B 2 Perror */
3642  pressure.PresTotlError*100.,
3643  /*C 3 Pcurrent */
3644  pressure.PresTotlCurr,
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 */
3648  pressure.PresTotlInit + pressure.PresInteg - pressure.pinzon,
3649  /*E 5 pgas (0) */
3650  pressure.PresTotlInit,
3651  /*F 6 Pgas */
3652  pressure.PresGasCurr,
3653  /*G 7 Pram */
3654  pressure.PresRamCurr,
3655  /*H 8 P rad in lines */
3656  pressure.pres_radiation_lines_curr,
3657  /*I 9 Pinteg subtract continuum rad pres which has already been added on */
3658  pressure.PresInteg - pressure.pinzon,
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 */
3662  timesc.sound_speed_adiabatic/1e5,
3663  /* the magnetic pressure */
3664  magnetic.pressure ,
3665  /* the local turbulent velocity in km/s */
3666  DoppVel.TurbVel/1e5 ,
3667  /* turbulent pressure */
3668  pressure.PresTurbCurr*DoppVel.lgTurb_pressure ,
3669  /* gravitational pressure */
3670  pressure.IntegRhoGravity,
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
3675  pressure.PresIntegElecThin-pressure.pinzon_PresIntegElecThin,
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 */
3686  radius.depth_mid_zone ,
3687  pressure.PresTotlCurr,
3688  pressure.pres_radiation_lines_curr/SDIV(pressure.PresTotlCurr) );
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",
3702  nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
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",
3713  nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
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 );
3730  cdEXIT(EXIT_FAILURE);
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 ,
3763  secondaries.csupra[ipHYDROGEN][0],
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)),
3780  rfield.ConEmitLocal[nzone][j]/rfield.widflx(j),
3781  opac.opacity_abs[j],
3782  rfield.ConSourceFcnLocal[nzone][j],
3783  rfield.ConSourceFcnLocal[nzone][j]/plankf(j),
3784  safe_div(rfield.ConSourceFcnLocal[nzone][j],rfield.flux[0][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",
3840  radius.depth_mid_zone,
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",
3892  nzone, phycon.TeInit, phycon.TeProp, phycon.te,
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",
3908  radius.Radius_mid_zone,
3909  radius.depth_mid_zone,
3910  wind.windv,
3911  wind.AccelTotalOutward,
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" );
3931  cdEXIT(EXIT_FAILURE);
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" );
3947  cdEXIT(EXIT_FAILURE);
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" );
3963  cdEXIT(EXIT_FAILURE);
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" );
3979  cdEXIT(EXIT_FAILURE);
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" );
3995  cdEXIT(EXIT_FAILURE);
3996  }
3997  }
3998  }
3999  else if( strcmp(save.chSave[ipPun],"XLNR") == 0 )
4000  {
4001  /* reflected lines */
4002  ASSERT( grid.lgOutputTypeOn[7] );
4003 
4004  if( lgLastOnly )
4005  {
4006  if( grid.lgGrid )
4007  saveFITSfile( save.params[ipPun].ipPnunit, 7 );
4008  else
4009  {
4010  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4011  cdEXIT(EXIT_FAILURE);
4012  }
4013  }
4014  }
4015  else if( strcmp(save.chSave[ipPun],"XLNO") == 0 )
4016  {
4017  /* outward lines */
4018  ASSERT( grid.lgOutputTypeOn[6] );
4019 
4020  if( lgLastOnly )
4021  {
4022  if( grid.lgGrid )
4023  saveFITSfile( save.params[ipPun].ipPnunit, 6 );
4024  else
4025  {
4026  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4027  cdEXIT(EXIT_FAILURE);
4028  }
4029  }
4030  }
4031  else if( strcmp(save.chSave[ipPun],"XREF") == 0 )
4032  {
4033  /* total reflected, lines and continuum */
4034  ASSERT( grid.lgOutputTypeOn[9] );
4035 
4036  if( lgLastOnly )
4037  {
4038  if( grid.lgGrid )
4039  saveFITSfile( save.params[ipPun].ipPnunit, 9 );
4040  else
4041  {
4042  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4043  cdEXIT(EXIT_FAILURE);
4044  }
4045  }
4046  }
4047  else if( strcmp(save.chSave[ipPun],"XTOT") == 0 )
4048  {
4049  /* total spectrum, reflected plus transmitted */
4050  ASSERT( grid.lgOutputTypeOn[0] );
4051 
4052  if( lgLastOnly )
4053  {
4054  if( grid.lgGrid )
4055  saveFITSfile( save.params[ipPun].ipPnunit, 0 );
4056  else
4057  {
4058  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4059  cdEXIT(EXIT_FAILURE);
4060  }
4061  }
4062  }
4063  else if( strcmp(save.chSave[ipPun],"XTRN") == 0 )
4064  {
4065  /* total outward, lines and continuum */
4066  ASSERT( grid.lgOutputTypeOn[8] );
4067 
4068  if( lgLastOnly )
4069  {
4070  if( grid.lgGrid )
4071  saveFITSfile( save.params[ipPun].ipPnunit, 8 );
4072  else
4073  {
4074  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4075  cdEXIT(EXIT_FAILURE);
4076  }
4077  }
4078  }
4079  else if( strcmp(save.chSave[ipPun],"XSPM") == 0 )
4080  {
4081  /* exp(-tau) to the illuminated face */
4082  ASSERT( grid.lgOutputTypeOn[10] );
4083 
4084  if( lgLastOnly )
4085  {
4086  if( grid.lgGrid )
4087  saveFITSfile( save.params[ipPun].ipPnunit, 10 );
4088  else
4089  {
4090  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4091  cdEXIT(EXIT_FAILURE);
4092  }
4093  }
4094  }
4095  // termination of second set of nested if's
4096  // error if we have not matched key
4097  /* there are a few "save" commands that are handled elsewhere
4098  * save dr is an example. These will have lgRealSave set false */
4099  // lgNoHitFirstBranch says did not find in previous nest of if's
4100  else if( save.lgRealSave[ipPun] && lgNoHitFirstBranch )
4101  {
4102  /* this is insanity, internal flag set in ParseSave not seen here */
4103  fprintf( ioQQQ, " PROBLEM DISASTER SaveDo does not recognize flag %4.4s set by ParseSave. This is impossible.\n",
4104  save.chSave[ipPun] );
4105  TotalInsanity();
4106  }
4107 
4108  /* print special hash string to separate out various iterations
4109  * chTime is LAST on last iteration
4110  * save.lgHashEndIter flag is true by default, set false
4111  * with "no hash" keyword on save command
4112  * save.lg_separate_iterations is true by default, set false
4113  * when save time dependent calcs since do not want special
4114  * character between time steps
4115  * grid.lgGrid is only true when doing a grid of calculations */
4116  if( lgLastOnly &&
4117  !(iterations.lgLastIt && !grid.lgGrid ) &&
4118  save.lgHashEndIter[ipPun] &&
4119  save.lg_separate_iterations[ipPun] &&
4120  !save.lgFITS[ipPun] )
4121  {
4122  if( dynamics.lgTimeDependentStatic && strcmp( save.chHashString , "TIME_DEP" )==0 )
4123  {
4124  fprintf( save.params[ipPun].ipPnunit, "\"time=%f\n",
4125  dynamics.time_elapsed );
4126  }
4127  else if( strcmp( save.chHashString , "\n" )==0 )
4128  {
4129  fprintf( save.params[ipPun].ipPnunit, "%s\n",
4130  save.chHashString );
4131  }
4132  else
4133  {
4134  fprintf( save.params[ipPun].ipPnunit, "%s",
4135  save.chHashString );
4136  if( grid.lgGrid && ( iterations.lgLastIt || lgAbort ) )
4137  fprintf( save.params[ipPun].ipPnunit, " GRID_DELIMIT -- grid%09ld",
4138  optimize.nOptimiz );
4139  fprintf( save.params[ipPun].ipPnunit, "\n" );
4140  }
4141  }
4142  if( save.lgFLUSH )
4143  fflush( save.params[ipPun].ipPnunit );
4144  }
4145  }
4146  return;
4147 }
4148 
4149 /*SaveLineIntensity produce the 'save lines intensity' output */
4150 STATIC void SaveLineIntensity(FILE * ioPUN, long int ipPun , realnum Threshold )
4151 {
4152  long int i;
4153 
4154  DEBUG_ENTRY( "SaveLineIntensity()" );
4155 
4156  /* used to save out all the emission line intensities
4157  * first initialize the line image reader */
4158 
4159  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4160  input.echo(ioPUN);
4161 
4162  /* now print any cautions or warnings */
4163  cdWarnings( ioPUN);
4164  cdCautions( ioPUN);
4165  fprintf( ioPUN, "zone=%5ld\n", nzone );
4166  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4167  fprintf( ioPUN, "begin emission lines\n" );
4168 
4169 
4170  // check whether intrinsic or emergent line emissivity
4171  bool lgEmergent = false;
4172  if( save.punarg[ipPun][0] > 0 )
4173  lgEmergent = true;
4174 
4175  /* only save non-zero intensities */
4176  fixit("value of lgEmergent isn't consistent with lgEmergent");
4177  SaveLineResults slr(ioPUN,save.lgEmergent[ipPun]);
4178 
4179  for( i=0; i < LineSave.nsum; i++ )
4180  {
4181  // Threshold is zero by default on save line intensity,
4182  // all option sets to negative number so that we report all lines
4183  if( LineSave.lines[i].SumLine(lgEmergent) > Threshold )
4184  {
4185  slr.save(&LineSave.lines[i]);
4186  }
4187  }
4188 
4189  slr.flush();
4190 
4191  fprintf( ioPUN, " \n" );
4192  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4193 
4194  return;
4195 }
4196 
4197 /* lgSaveOpticalDepths true says save optical depths */
4198 static bool lgPopsFirstCall , lgSaveOpticalDepths;
4199 
4200 /*SaveLineStuff save optical depths or source functions for all transferred lines */
4201 STATIC void SaveLineStuff(
4202  FILE * ioPUN,
4203  const char *chJob ,
4204  realnum xLimit )
4205 {
4206  DEBUG_ENTRY( "SaveLineStuff()" );
4207 
4208  /*find out which job this is and set a flag to use later */
4209  if( strcmp( &*chJob , "optical" ) == 0 )
4210  {
4211  /* save line optical depths */
4212  lgSaveOpticalDepths = true;
4213  lgPopsFirstCall = false;
4214  }
4215  else if( strcmp( &*chJob , "populat" ) == 0 )
4216  {
4217  static bool lgFirst=true;
4218  lgSaveOpticalDepths = false;
4219  /* level population information */
4220  if( lgFirst )
4221  {
4222  lgPopsFirstCall = true;
4223  fprintf(ioPUN,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n");
4224  lgFirst = false;
4225  }
4226  else
4227  {
4228  lgPopsFirstCall = false;
4229  }
4230  }
4231  else
4232  {
4233  fprintf( ioQQQ, " insane job in SaveLineStuff =%s\n",
4234  &*chJob );
4235  cdEXIT(EXIT_FAILURE);
4236  }
4237 
4238  long index = 0;
4239  /* loop over all lines, calling put1Line to create info (routine located below) */
4240  /* hydrogen like lines */
4241  /* >>chng 02 may 16, had been explicit H and He-like loops */
4242  for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
4243  {
4244  for( long nelem=ipISO; nelem < LIMELM; nelem++ )
4245  {
4246  if( dense.lgElmtOn[nelem] )
4247  {
4248  /* 06 aug 28, from numLevels_max to _local. */
4249  for( long ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
4250  {
4251  for( long ipLo=0; ipLo <ipHi; ipLo++ )
4252  {
4253  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
4254  continue;
4255 
4256  ++index;
4257  Save1Line( iso_sp[ipISO][nelem].trans(ipHi,ipLo), ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) );
4258  }
4259  }
4260  /* also do extra Lyman lines if optical depths are to be done,
4261  * these are line that are included only for absorption, not in the
4262  * model atoms */
4263  if( lgSaveOpticalDepths )
4264  {
4265  /* >>chng 02 aug 23, for he-like, had starting on much too high a level since
4266  * index was number of levels - caught by Adrian Turner */
4267  /* now output extra line lines, starting one above those already done above */
4268  /*for( ipHi=iso_sp[ipISO][nelem].numLevels_max; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )*/
4269  /* 06 aug 28, from numLevels_max to _local. */
4270  for( long ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
4271  {
4272  ++index;
4273  Save1Line( ExtraLymanLines[ipISO][nelem][ipExtraLymanLines[ipISO][nelem][ipHi]], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) );
4274  }
4275  }
4276  }
4277  }
4278  }
4279 
4280  for( long i=0; i < nWindLine; i++ )
4281  {
4282  if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
4283  {
4284  ++index;
4285  Save1Line( TauLine2[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
4286  }
4287  }
4288 
4289  for( size_t i=0; i < UTALines.size(); i++ )
4290  {
4291  ++index;
4292  Save1Line( UTALines[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) );
4293  }
4294 
4295  /* save optical depths of H2 lines */
4296  h2.H2_PunchLineStuff( ioPUN , xLimit , index);
4297 
4298  /* data base lines */
4299  for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
4300  {
4301  realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
4302  for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
4303  em != dBaseTrans[ipSpecies].Emis().end(); ++em)
4304  {
4305  ++index;
4306  Save1Line( (*em).Tran(), ioPUN, xLimit, index, DopplerWidth );
4307  }
4308  }
4309 
4310  /*fprintf(ioPUN, "##################################\n"); */
4311  fprintf( ioPUN , "%s\n",save.chHashString );
4312  return;
4313 }
4314 
4315 /*Save1Line called by SaveLineStuff to produce output for one line */
4316 void Save1Line( const TransitionProxy& t , FILE * ioPUN , realnum xLimit , long index, realnum DopplerWidth )
4317 {
4318 
4319  if( lgSaveOpticalDepths )
4320  {
4321  /* optical depths, no special first time, only print them */
4322  if( t.Emis().TauIn() >= xLimit )
4323  {
4324  /* label like "C 4" or "H 1" */
4325  fprintf( ioPUN, "%-*.*s\t",CHARS_SPECIES, CHARS_SPECIES, chIonLbl(t).c_str());
4326 
4327  /* print wavelengths, either line in main printout labels,
4328  * or in various units in exponential notation - prt_wl is in prt.c */
4329  if( strcmp( save.chConSavEnr[save.ipConPun], "labl" )== 0 )
4330  {
4331  prt_wl( ioPUN , t.WLAng() );
4332  }
4333  else
4334  {
4335  /* this converts energy in Rydbergs into any of the other units */
4336  fprintf( ioPUN , "%.7e", AnuUnit((realnum)(t.EnergyRyd())) );
4337  }
4338  /* print the optical depth */
4339  fprintf( ioPUN , "\t%.3f", t.Emis().TauIn() );
4340  /* damping constant */
4341  fprintf(ioPUN, "\t%.3e",
4342  t.Emis().dampXvel() / DopplerWidth );
4343  fprintf(ioPUN, "\n");
4344  }
4345  }
4346  else if( lgPopsFirstCall )
4347  {
4348  /* first call to line populations, print atomic parameters and indices */
4349  fprintf(ioPUN, "%li\t%s" , index , chLineLbl(t).c_str());
4350  /* stat weights */
4351  fprintf(ioPUN, "\t%.0f\t%.0f",
4352  (*t.Lo()).g() ,(*t.Hi()).g());
4353  /* energy difference, gf */
4354  fprintf(ioPUN, "\t%.2f\t%.3e",
4355  t.EnergyWN() ,t.Emis().gf());
4356  fprintf(ioPUN, "\n");
4357  }
4358  else
4359  {
4360  /* not first call, so do level populations and indices defined above */
4361  if( (*t.Hi()).Pop() > xLimit )
4362  {
4363  /* >>chng 05 may 08, add abundances, which for iso-seq species is
4364  * the density of the parent ion, for other lines, is unity.
4365  * had not been included so pops for iso seq were rel to parent ion.
4366  * caught by John Everett */
4367  /* multiplication by abundance no longer necessary since iso pops now denormalized */
4368  fprintf(ioPUN,"%li\t%.2e\t%.2e\n", index, (*t.Lo()).Pop(), (*t.Hi()).Pop() );
4369  }
4370  }
4371 }
4372 
4373 /*SaveNewContinuum produce the 'save new continuum' output */
4374 STATIC void SaveNewContinuum(FILE * ioPUN )
4375 {
4376  DEBUG_ENTRY( "SaveNewContinuum()" );
4377  long int ipLo, ipHi,
4378  j ,
4379  ncells;
4380 
4381  double wllo, wlhi;
4382 
4383  double *cont_incid,
4384  *cont_atten,
4385  *diffuse_in,
4386  *diffuse_out,
4387  *emis_lines_out,
4388  *emis_lines_in;
4389 
4390  /* get the low limit */
4391  wllo = rfield.anu(0);
4392  /* get high-energy limit */
4393  wlhi = rfield.anu(rfield.nflux-1);
4394  /* use native continuum mesh */
4395  ipLo = ipoint(wllo)-1;
4396  ipHi = ipoint(wlhi)-1;
4397  ncells = ipHi - ipLo + 1;
4398 
4399  /* now allocate the space */
4400  cont_incid = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4401  cont_atten = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4402  diffuse_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4403  diffuse_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4404  emis_lines_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4405  emis_lines_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4406  /*emis_lines_pump_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));
4407  emis_lines_pump_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));*/
4408 
4409  fixit("all of these should be rerouted to cdSPEC2, but units not right");
4410  /* all of these should be rerouted to cdSPEC2, but units not right
4411  * need ergs multiplied for one, and continuum and lines may not be added correctly.
4412  * Goal is to abandon current cdSPEC and replace it with current cdSPEC2 */
4413 
4414 #if 1
4415  /* for cdSPEC the energy vector is the lower edge of the energy cell */
4416  /* get incident continuum */
4417  cdSPEC( 1 , ncells , cont_incid );
4418  /* get attenuated incident continuum */
4419  cdSPEC( 2 , ncells , cont_atten );
4420  /* get diffuse continuous emission, reflected */
4421  cdSPEC( 5 , ncells , diffuse_in );
4422  /* get continuous emission outward direction */
4423  cdSPEC( 4 , ncells , diffuse_out );
4425  /* get all outward lines */
4426  cdSPEC( 6 , ncells , emis_lines_out );
4427  /* get all reflected lines */
4428  cdSPEC( 7 , ncells , emis_lines_in );
4429 #else
4430  cdSPEC2( 1, rfield.nflux, 0, rfield.nflux - 1, cont_incid );
4431  /* get attenuated incident continuum */
4432  cdSPEC2( 2, rfield.nflux, 0, rfield.nflux - 1, cont_atten );
4433  /* get diffuse continuous emission, reflected */
4434  cdSPEC2( 5, rfield.nflux, 0, rfield.nflux - 1, diffuse_in );
4435  /* get continuous emission outward direction */
4436  cdSPEC2( 4, rfield.nflux, 0, rfield.nflux - 1, diffuse_out );
4437  /* get all outward lines */
4438  cdSPEC2( 6, rfield.nflux, 0, rfield.nflux - 1, emis_lines_out );
4439  /* get all reflected lines */
4440  cdSPEC2( 7, rfield.nflux, 0, rfield.nflux - 1, emis_lines_in );
4441 #endif
4442 
4443  /* for this example we will do a wavelength range */
4444  for( j=0; j<ncells-1; ++j )
4445  {
4446  /* photon energy in appropriate energy or wavelength units */
4447  fprintf( ioPUN,"%.5e\t", AnuUnit((realnum)(rfield.anu(j+ipLo)) ) );
4448  fprintf( ioPUN,"%.3e\t", cont_incid[j] );
4449  fprintf( ioPUN,"%.3e\t", cont_atten[j] );
4450  fprintf( ioPUN,"%.3e\t", diffuse_in[j]+diffuse_out[j] );
4451  fprintf( ioPUN,"%.3e",
4452  emis_lines_out[j]+emis_lines_in[j]/*+emis_lines_pump_out[j]+emis_lines_pump_in[j]*/ );
4453  fprintf( ioPUN,"\n" );
4454  }
4455 
4456  free(cont_incid);
4457  free(diffuse_in);
4458  free(diffuse_out);
4459  free(cont_atten);
4460  free(emis_lines_out);
4461  free(emis_lines_in);
4462  /*free(emis_lines_pump_out);
4463  free(emis_lines_pump_in );*/
4464 }
4465 
4466 /* save AGN3 hemiss, for Chapter 4, routine is below */
4467 STATIC void AGN_Hemis(FILE *ioPUN )
4468 {
4469  const int NTE = 4;
4470  realnum te[NTE] = {5000., 10000., 15000., 20000.};
4471  realnum *agn_continuum[NTE];
4472  double TempSave = phycon.te;
4473  long i , j;
4474 
4475  DEBUG_ENTRY( "AGN_Hemis()" );
4476 
4477  /* make table of continuous emission at various temperatuers */
4478  /* first allocate space */
4479  for( i=0;i<NTE; ++i)
4480  {
4481  agn_continuum[i] = (realnum *)MALLOC((unsigned)rfield.nflux*sizeof(realnum) );
4482 
4483  /* set the next temperature */
4484  /* recompute everything at this new temp */
4485  TempChange(te[i] , true);
4486  /* converge the pressure-temperature-ionization solution for this zone */
4487  ConvPresTempEdenIoniz();
4488 
4489  /* now get the thermal emission */
4490  RT_diffuse();
4491  for(j=0;j<rfield.nflux; ++j )
4492  {
4493  agn_continuum[i][j] = rfield.ConEmitLocal[nzone][j]/(realnum)dense.eden/
4494  (dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHELIUM][1] + dense.xIonDense[ipHELIUM][2] );
4495  }
4496  }
4497 
4498  /* print title for line */
4499  fprintf(ioPUN,"wl");
4500  for( i=0;i<NTE; ++i)
4501  {
4502  fprintf(ioPUN,"\tT=%.0f",te[i]);
4503  }
4504  fprintf( ioPUN , "\tcont\n");
4505 
4506  /* not print all n temperatures across a line */
4507  for(j=0;j<rfield.nflux; ++j )
4508  {
4509  fprintf( ioPUN , "%12.5e",
4510  AnuUnit(rfield.anu(j)) );
4511  /* loop over the temperatures, and for each, calculate a continuum */
4512  for( i=0;i<NTE; ++i)
4513  {
4514  fprintf(ioPUN,"\t%.3e",agn_continuum[i][j]*rfield.anu2(j)*EN1RYD/rfield.widflx(j));
4515  }
4516  /* cont label and end of line*/
4517  fprintf( ioPUN , "\t%s\n" , rfield.chContLabel[j].c_str());
4518  }
4519 
4520  /* now free the continua */
4521  for( i=0;i<NTE; ++i)
4522  {
4523  free( agn_continuum[i] );
4524  }
4525 
4526  /* Restore temperature stored before this routine was called */
4527  /* and force update */
4528  TempChange(TempSave , true);
4529 
4530  fprintf( ioQQQ, "AGN_Hemis - result of save AGN3 hemis - I have left the code in a disturbed state, and will now exit.\n");
4531  cdEXIT(EXIT_FAILURE);
4532 }
4533 
4534 /*SaveResults save results from save results command */
4535 /*SaveResults1Line do single line of output for the save results and save line intensity commands */
4536 STATIC void SaveResults(FILE* ioPUN)
4537 {
4538  long int i , nelem , ion;
4539 
4540  DEBUG_ENTRY( "SaveResults()" );
4541 
4542  /* used to save out line intensities, optical depths,
4543  * and column densities */
4544 
4545  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4546  input.echo(ioPUN);
4547 
4548  /* first print any cautions or warnings */
4549  cdWarnings(ioPUN);
4550  cdCautions(ioPUN);
4551  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4552 
4553  fprintf( ioPUN, "C*OPTICAL DEPTHS ELECTRON=%10.3e\n", opac.telec );
4554 
4555  fprintf( ioPUN, "BEGIN EMISSION LINES\n" );
4556  SaveLineResults slr(ioPUN,0);
4557 
4558  for( i=0; i < LineSave.nsum; i++ )
4559  {
4560  if( LineSave.lines[i].SumLine(0) > 0. )
4561  {
4562  slr.save(&LineSave.lines[i]);
4563  }
4564  }
4565 
4566  slr.flush();
4567 
4568  fprintf( ioPUN, " \n" );
4569 
4570  fprintf( ioPUN, "BEGIN COLUMN DENSITIES\n" );
4571 
4572  /* this dumps out the whole array,*/
4573  /* following loop relies on LIMELM being 30, assert it here in case
4574  * this is ever changed */
4575  ASSERT( LIMELM == 30 );
4576  /* this order of indices is to keep 30 as the fastest variable,
4577  * and the 32 (LIMELM+1) as the slower one */
4578  for( nelem=0; nelem<LIMELM; nelem++ )
4579  {
4580  for(ion=0; ion < nelem+1; ion++)
4581  {
4582  fprintf( ioPUN, " %10.3e", mean.xIonMean[0][nelem][ion][0] );
4583  /* throw line feed every 10 numbers */
4584  if( nelem==9|| nelem==19 || nelem==29 )
4585  {
4586  fprintf( ioPUN, "\n" );
4587  }
4588  }
4589  }
4590 
4591  fprintf( ioPUN, "END OF RESULTS\n" );
4592  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4593  return;
4594 }
4595 
4596 namespace
4597 {
4598  /*SaveResults1Line do single line of output for the save results and save line intensity commands */
4599  /* the number of emission lines across one line of printout */
4600  void SaveLineResults::save(
4601  /* 4 char + null string */
4602  const LinSv *line)
4603  {
4604 
4605  DEBUG_ENTRY( "SaveLineResults::save()" );
4606 
4607  /* if LineWidth is changed then change format in write too */
4608 
4609  /* save results in array so that they can be printed when done */
4610  m_lines[ipLine] = line;
4611 
4612  /* now increment the counter and then check if we have filled the line,
4613  * and so should print it */
4614  ++ipLine;
4615  /* do print now if we are in column mode (one line per line) or if we have filled up
4616  * the line */
4617  if( ( strcmp(::save.chPunRltType,"column") == 0 ) || ipLine == LINEWIDTH )
4618  {
4619  /* "array " is usual array 6 wide, "column" is one line per line */
4620  flush();
4621  }
4622  }
4623 }
4624 
4625 /*SaveGaunts called by save gaunts command to output Gaunt factors */
4626 STATIC void SaveGaunts(FILE* ioPUN)
4627 {
4628  DEBUG_ENTRY( "SaveGaunts()" );
4629 
4630  static const int NENR_GAUNT = 33;
4631  static const int NTE_GAUNT = 20;
4632 
4633  double ener[NENR_GAUNT], ste[NTE_GAUNT], g[NENR_GAUNT][NTE_GAUNT];
4634 
4635  /* this routine is called from the PUNCH GAUNTS command
4636  * it drives the Gaunt factor routine to save gaunts over full range */
4637 
4638  for( int i=0; i < NTE_GAUNT; i++ )
4639  ste[i] = 0.5f*(i+1);
4640 
4641  for( int i=0; i < NENR_GAUNT; i++ )
4642  ener[i] = 0.5f*i - 9.f;
4643 
4644  for( int charge=1; charge <= LIMELM; charge++ )
4645  {
4646  /* energy is log of energy */
4647  for( int ite=0; ite < NTE_GAUNT; ite++ )
4648  {
4649  for( int j=0; j < NENR_GAUNT; j++ )
4650  {
4651  g[j][ite] = t_gaunt::Inst().gauntff( charge, exp10(ste[ite]), exp10(ener[j]) );
4652  }
4653  }
4654 
4655  /* now save out the results */
4656  fprintf( ioPUN, "\tlg(nu)\\lg(Te)" );
4657  for( int i=1; i <= NTE_GAUNT; i++ )
4658  fprintf( ioPUN, "\t%.3e", ste[i-1] );
4659  fprintf( ioPUN, "\n" );
4660 
4661  for( int j=0; j < NENR_GAUNT; j++ )
4662  {
4663  fprintf( ioPUN, "\t%10.3e", ener[j] );
4664  for( int ite=0; ite < NTE_GAUNT; ite++ )
4665  fprintf( ioPUN, "\t%.3e", g[j][ite] );
4666  fprintf( ioPUN, "\n" );
4667  }
4668 
4669  fprintf( ioPUN, "\tlg(nu)/lg(Te)" );
4670  for( int i=0; i < NTE_GAUNT; i++ )
4671  fprintf( ioPUN, "\t%.3e", ste[i] );
4672  fprintf( ioPUN, "\n\n" );
4673 
4674  fprintf( ioPUN, "Below is log(gamma^2), log(u), gff\n" );
4675  /* print log(gamma2), log(u) instead of temp and energy. */
4676 
4677  double z = log10((double)charge);
4678 
4679  for( int i=0; i < NTE_GAUNT; i++ )
4680  {
4681  for( int j=0; j < NENR_GAUNT; j++ )
4682  {
4683  fprintf( ioPUN, "\t%10.3e\t%10.3e\t%10.3e\n",
4684  2.*z + log10(TE1RYD) - ste[i],
4685  log10(TE1RYD) + ener[j] - ste[i],
4686  g[j][i] );
4687  }
4688  }
4689  fprintf( ioPUN, "end of charge = %i\n", charge );
4690  fprintf( ioPUN, "****************************\n" );
4691  }
4692 }
4693 
4694 void SaveGrid(FILE* pnunit, exit_type status)
4695 {
4696  DEBUG_ENTRY( "SaveGrid()" );
4697 
4698  if( pnunit == NULL )
4699  return;
4700 
4701  if( optimize.nOptimiz == 0 )
4702  {
4703  /* start of line gives abort and warning summary */
4704  fprintf( pnunit, "#Index\tFailure?\tWarnings?\tExit code\t#rank\t#seq" );
4705  /* print start of each variable command line */
4706  for( int i=0; i < grid.nintparm; i++ )
4707  {
4708  char chStr[10];
4709  strncpy( chStr, optimize.chVarFmt[i], 9 );
4710  /* make sure this small bit of string is terminated */
4711  chStr[9] = '\0';
4712  fprintf( pnunit, "\t%s", chStr );
4713  }
4714  fprintf( pnunit, "\tgrid parameter string\n" );
4715  }
4716  /* abort / warning summary for this sim */
4717  bool lgNoFailure = ( status == ES_SUCCESS || status == ES_WARNINGS );
4718  fprintf( pnunit, "%9.9ld\t%c\t%c\t%20s\t%ld\t%ld",
4719  optimize.nOptimiz,
4720  TorF(!lgNoFailure),
4721  TorF(warnings.lgWarngs),
4722  cpu.i().chExitStatus(status).c_str(),
4723  cpu.i().nRANK(),
4724  grid.seqNum );
4725  /* the grid parameters */
4726  char chGridParam[INPUT_LINE_LENGTH];
4727  char chStringHold[100];
4728  sprintf( chStringHold, "%f", grid.interpParameters[optimize.nOptimiz][0] );
4729  strcpy( chGridParam, chStringHold );
4730  for( int j=0; j < grid.nintparm; j++ )
4731  {
4732  if( j > 0 )
4733  {
4734  sprintf( chStringHold, ", %f", grid.interpParameters[optimize.nOptimiz][j] );
4735  strcat( chGridParam, chStringHold );
4736  }
4737  fprintf( pnunit, "\t%f", grid.interpParameters[optimize.nOptimiz][j] );
4738  }
4739  fprintf( pnunit, "\t%s\n", chGridParam );
4740 }
SaveLineStuff
STATIC void SaveLineStuff(FILE *ioPUN, const char *chJob, realnum xLimit)
Definition: save_do.cpp:4201
t_cpu_i::nRANK
long nRANK() const
Definition: cpu.h:392
AGN_He1_CS
void AGN_He1_CS(FILE *ioPun)
Definition: iso_solve.cpp:503
monitor_results.h
t_save::lgPunLstIter
bool lgPunLstIter[LIMPUN]
Definition: save.h:374
t_conv::getCounterZone
long getCounterZone(const long type) const
Definition: conv.h:340
t_secondaries::x12tot
realnum x12tot
Definition: secondaries.h:65
t_pressure::PresTotlInit
double PresTotlInit
Definition: pressure.h:52
t_mesh::emm
double emm() const
Definition: mesh.h:84
t_rfield::ConSourceFcnLocal
realnum ** ConSourceFcnLocal
Definition: rfield.h:144
t_ionbal::DR_Badnell_rate_coef
double ** DR_Badnell_rate_coef
Definition: ionbal.h:197
HydroRecCool
double HydroRecCool(long int n, long int ipZ)
Definition: hydroreccool.cpp:10
mole_global
t_mole_global mole_global
Definition: mole.cpp:7
t_rfield::fine_opt_depth
realnum * fine_opt_depth
Definition: rfield.h:391
t_rfield::fine_anu
realnum * fine_anu
Definition: rfield.h:393
t_conv::hist_pres_density
vector< double > hist_pres_density
Definition: conv.h:293
TexcLine
double TexcLine(const TransitionProxy &t)
Definition: transition.cpp:204
t_hmi::H2_Solomon_dissoc_rate_used_H2g
double H2_Solomon_dissoc_rate_used_H2g
Definition: hmi.h:103
t_radius::Radius
double Radius
Definition: radius.h:31
t_radius::depth
double depth
Definition: radius.h:31
t_save::punarg
realnum punarg[LIMPUN][3]
Definition: save.h:357
GrainVar::dstab
vector< double > dstab
Definition: grainvar.h:527
t_thermal::htot
double htot
Definition: thermal.h:169
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:44
t_rfield::line_count
long int * line_count
Definition: rfield.h:61
t_opac::ipElement
long int ipElement[LIMELM][LIMELM][7][3]
Definition: opacity.h:222
Wind::dvdr
realnum dvdr
Definition: wind.h:21
Wind::AccelLine
realnum AccelLine
Definition: wind.h:61
thermal
t_thermal thermal
Definition: thermal.cpp:6
prt_LineLabels
void prt_LineLabels(FILE *ioOUT, bool lgPrintAll)
Definition: prt.cpp:112
SaveSpecial
void SaveSpecial(FILE *io, const char *chTime)
Definition: save_special.cpp:12
t_LineSave::SortWL
vector< size_t > SortWL
Definition: lines.h:134
SaveHeat
void SaveHeat(FILE *io)
Definition: heat_save.cpp:22
chIonLbl
string chIonLbl(const TransitionProxy &t)
Definition: transition.cpp:230
ipMAGNESIUM
const int ipMAGNESIUM
Definition: cddefines.h:359
CoolSave
void CoolSave(FILE *io, const char chJob[])
Definition: cool_save.cpp:20
TransitionList::size
size_t size(void) const
Definition: transition.h:331
colden
t_colden colden
Definition: colden.cpp:5
is_odd
bool is_odd(int j)
Definition: cddefines.h:757
t_opac::opacity_abs
double * opacity_abs
Definition: opacity.h:103
t_pressure::PresInteg
realnum PresInteg
Definition: pressure.h:69
ipPun
STATIC long int ipPun
Definition: save_do.cpp:721
null_mole
molecule * null_mole
Definition: mole_species.cpp:50
t_iso_sp::st
qList st
Definition: iso.h:482
t_dynamics::Cool
double Cool()
Definition: dynamics.cpp:2205
t_phycon::te03
double te03
Definition: phycon.h:58
UTALines
TransitionList UTALines("UTALines",&AnonStates)
t_opac::albedo
double * albedo
Definition: opacity.h:112
exp10
double exp10(double x)
Definition: cddefines.h:1383
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:65
t_conv::hist_pres_error
vector< double > hist_pres_error
Definition: conv.h:293
chLineLbl
string chLineLbl(const TransitionProxy &t)
Definition: transition.h:599
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:1067
t_mesh::widflx
double widflx(size_t i) const
Definition: mesh.h:147
t_dense::EdenHCorr
double EdenHCorr
Definition: dense.h:227
input
t_input input
Definition: input.cpp:12
t_radius::PI4_rinner_sq
double PI4_rinner_sq
Definition: radius.h:31
t_grid::lgGrid
bool lgGrid
Definition: grid.h:42
opac
t_opac opac
Definition: opacity.cpp:5
t_save::ResolutionAbs
realnum ResolutionAbs
Definition: save.h:478
t_mole_global::num_calc
int num_calc
Definition: mole.h:362
t_rfield::flux
realnum ** flux
Definition: rfield.h:70
struc
t_struc struc
Definition: struc.cpp:6
t_rfield::DiffuseLineEmission
realnum * DiffuseLineEmission
Definition: rfield.h:195
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
t_iso_sp::RadRec_caseB
double RadRec_caseB
Definition: iso.h:544
hyperfine
t_hyperfine hyperfine
Definition: hyperfine.cpp:5
t_ADfA::ph1
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat_adfa.h:61
t_rfield::chContLabel
vector< string > chContLabel
Definition: rfield.h:215
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:246
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
t_cpu::i
t_cpu_i & i()
Definition: cpu.h:415
NISO
const int NISO
Definition: cddefines.h:310
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:49
t_hcmap::MapZone
long int MapZone
Definition: hcmap.h:20
t_mesh::mesh_md5sum
string mesh_md5sum() const
Definition: mesh.h:103
t_phycon::EnergyIonization
double EnergyIonization
Definition: phycon.h:41
GrainVar::GasCoolColl
double GasCoolColl
Definition: grainvar.h:546
t_rfield::outlin_noplot
realnum * outlin_noplot
Definition: rfield.h:191
save_line
void save_line(FILE *ip, const char *chDo, bool lgEmergent, long ipPun)
Definition: save_line.cpp:110
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:48
t_pressure::PresIntegElecThin
realnum PresIntegElecThin
Definition: pressure.h:75
t_pressure::PresRamCurr
double PresRamCurr
Definition: pressure.h:39
VERSION_TRNCON
static const long VERSION_TRNCON
Definition: save.h:16
cdCautions
void cdCautions(FILE *ioOUT)
Definition: cddrive.cpp:220
t_rfield::flux_total_incident
realnum ** flux_total_incident
Definition: rfield.h:201
TorF
char TorF(bool l)
Definition: cddefines.h:753
DynaSave
void DynaSave(FILE *ipPnunit, char chJob)
Definition: dynamics.cpp:2166
t_pressure::IntegRhoGravity
double IntegRhoGravity
Definition: pressure.h:83
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:46
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:355
t_conv::lgConvPres
bool lgConvPres
Definition: conv.h:194
magnetic
t_magnetic magnetic
Definition: magnetic.cpp:17
t_rfield::ipMaxBolt
long int ipMaxBolt
Definition: rfield.h:232
t_dynamics::lgTimeDependentStatic
bool lgTimeDependentStatic
Definition: dynamics.h:102
EmissionProxy::TauTot
realnum & TauTot() const
Definition: emission.h:478
t_iso_sp::nCollapsed_max
long int nCollapsed_max
Definition: iso.h:518
warnings
t_warnings warnings
Definition: warnings.cpp:11
cdWarnings
void cdWarnings(FILE *ioPNT)
Definition: cddrive.cpp:192
CHIANTI_Upsilon
double CHIANTI_Upsilon(long, long, long, long, double)
Definition: species2.cpp:900
conv
t_conv conv
Definition: conv.cpp:5
SaveGaunts
STATIC void SaveGaunts(FILE *ioPUN)
Definition: save_do.cpp:4626
t_optimize::nOptimiz
long int nOptimiz
Definition: optimize.h:250
HFLines
TransitionList HFLines("HFLines",&AnonStates)
t_mole_local::findrk
double findrk(const char buf[]) const
Definition: mole_reactions.cpp:4137
phycon
t_phycon phycon
Definition: phycon.cpp:6
t_phycon::EnthalpyDensity
double EnthalpyDensity
Definition: phycon.h:50
LineSave
t_LineSave LineSave
Definition: lines.cpp:9
t_DoppVel::TurbVel
realnum TurbVel
Definition: doppvel.h:21
SaveSpeciesBands
void SaveSpeciesBands(const long ipPun, const string &speciesLabel, const string &fileBands)
Definition: species_pseudo_cont.cpp:921
t_save::chHashString
char chHashString[INPUT_LINE_LENGTH]
Definition: save.h:398
ipHe2p1P
const int ipHe2p1P
Definition: iso.h:51
SaveLineData
NORETURN void SaveLineData(FILE *io)
Definition: save_linedata.cpp:20
matchGeneric
vector< genericState > matchGeneric(const char *chLabel, bool lgValidate)
Definition: generic_state.cpp:224
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:1095
ipRecNetEsc
const int ipRecNetEsc
Definition: cddefines.h:330
ion_recombAGN
void ion_recombAGN(FILE *io)
Definition: ion_recomb.cpp:216
ipFineCont
long ipFineCont(double energy_ryd)
Definition: cont_ipoint.cpp:165
CoolHeavy
t_CoolHeavy CoolHeavy
Definition: coolheavy.cpp:5
t_cpu_i::big_endian
bool big_endian() const
Definition: cpu.h:351
t_optimize::chVarFmt
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:267
t_geometry::covgeo
realnum covgeo
Definition: geometry.h:45
t_save::lgFLUSH
bool lgFLUSH
Definition: save.h:401
t_magnetic::pressure
double pressure
Definition: magnetic.h:41
t_rfield::outlin
realnum ** outlin
Definition: rfield.h:191
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
t_save::chPunRltType
char chPunRltType[7]
Definition: save.h:426
t_opac::opacity_sct
double * opacity_sct
Definition: opacity.h:106
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:829
PrtMeanIon
void PrtMeanIon(char chType, bool lgDensity, FILE *)
Definition: prt_meanion.cpp:11
t_geometry::FillFac
realnum FillFac
Definition: geometry.h:29
nzone
long int nzone
Definition: cddefines.cpp:14
t_hmi::HeatH2Dexc_used
double HeatH2Dexc_used
Definition: hmi.h:140
SaveParams::ipPnunit
FILE * ipPnunit
Definition: save.h:188
t_opac::TauScatFace
realnum * TauScatFace
Definition: opacity.h:99
t_called::lgTalk
bool lgTalk
Definition: called.h:12
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:43
DoppVel
t_DoppVel DoppVel
Definition: doppvel.cpp:5
dynamics
t_dynamics dynamics
Definition: dynamics.cpp:42
t_CoolHeavy::eebrm
double eebrm
Definition: coolheavy.h:18
t_save::lgRealSave
bool lgRealSave[LIMPUN]
Definition: save.h:288
t_phycon::TeProp
double TeProp
Definition: phycon.h:99
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:481
TauLine2
TransitionList TauLine2("TauLine2",&AnonStates)
molecule
Definition: mole.h:142
exit_type
exit_type
Definition: cddefines.h:142
PrtLinePres
void PrtLinePres(FILE *ioPRESSURE)
Definition: prt_linepres.cpp:16
MIN2
#define MIN2(a, b)
Definition: cddefines.h:807
t_phycon::te_eV
double te_eV
Definition: phycon.h:24
cdSPEC2
void cdSPEC2(int Option, long int nEnergy, long int ipLoEnergy, long int ipHiEnergy, realnum ReturnedSpectrum[])
t_pressure::PresTotlCurr
double PresTotlCurr
Definition: pressure.h:46
t_mesh::anu
double anu(size_t i) const
Definition: mesh.h:111
nSpecies
long int nSpecies
Definition: taulines.cpp:22
TempChange
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:31
t_save::lgEmergent
bool lgEmergent[LIMPUN]
Definition: save.h:291
t_LineSave::lines
vector< LinSv > lines
Definition: lines.h:132
ipSULPHUR
const int ipSULPHUR
Definition: cddefines.h:363
t_dynamics::lgStatic_completed
bool lgStatic_completed
Definition: dynamics.h:111
t_save::nsave
long int nsave
Definition: save.h:303
mole_save
void mole_save(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, bool lgCoef, double depth)
Definition: save_species.cpp:504
t_pressure::pres_radiation_lines_curr
double pres_radiation_lines_curr
Definition: pressure.h:61
GrainVar::GraphiteEmission
vector< realnum > GraphiteEmission
Definition: grainvar.h:581
dense
t_dense dense
Definition: global.cpp:15
Singleton< t_yield >::Inst
static t_yield & Inst()
Definition: cddefines.h:209
t_rfield::comup
double * comup
Definition: rfield.h:238
SaveDo
void SaveDo(const char *chTime)
Definition: save_do.cpp:817
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
t_rfield::nflux_with_check
long int nflux_with_check
Definition: rfield.h:51
GrainVar::dstsc
vector< double > dstsc
Definition: grainvar.h:528
atmdat_gaunt.h
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:391
t_dynamics::Heat
double Heat()
Definition: dynamics.cpp:2191
lgSaveOpticalDepths
static bool lgSaveOpticalDepths
Definition: save_do.cpp:4198
RT_diffuse
void RT_diffuse(void)
Definition: rt_diffuse.cpp:35
t_save::lgCumulative
bool lgCumulative[LIMPUN]
Definition: save.h:294
t_rfield::getCoarseTransCoef
const realnum * getCoarseTransCoef()
Definition: rfield.cpp:56
t_radius::PI4_Radius_sq
double PI4_Radius_sq
Definition: radius.h:31
t_ionbal::RateRecomTot
double ** RateRecomTot
Definition: ionbal.h:191
t_grid::interpParameters
realnum ** interpParameters
Definition: grid.h:32
t_magnetic::EnthalpyDensity
double EnthalpyDensity
Definition: magnetic.h:38
t_save::optname
string optname[LIMPUN]
Definition: save.h:360
wind
Wind wind
Definition: wind.cpp:5
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
t_geometry::lgSphere
bool lgSphere
Definition: geometry.h:34
t_conv::hist_temp_cool
vector< double > hist_temp_cool
Definition: conv.h:298
t_grid::nintparm
long nintparm
Definition: grid.h:57
iteration
long int iteration
Definition: cddefines.cpp:16
t_save::SpeciesBandFile
string SpeciesBandFile[LIMPUN]
Definition: save.h:495
t_rfield::otslin
realnum * otslin
Definition: rfield.h:185
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:858
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:556
t_save::LinEvery
long int LinEvery
Definition: save.h:462
SaveSpeciesOptDep
void SaveSpeciesOptDep(const long int ipPun, const string &speciesLabel)
Definition: save_species.cpp:418
t_radius::drad
double drad
Definition: radius.h:31
save_opacity
void save_opacity(FILE *io, long int np)
Definition: save_opacity.cpp:23
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:8
abund
t_abund abund
Definition: abund.cpp:5
GrainVar::TotalEden
double TotalEden
Definition: grainvar.h:531
Wind::AccelGravity
realnum AccelGravity
Definition: wind.h:49
t_save::lg_separate_iterations
bool lg_separate_iterations[LIMPUN]
Definition: save.h:319
t_rfield::ConEmitLocal
realnum ** ConEmitLocal
Definition: rfield.h:141
t_pressure::pinzon
realnum pinzon
Definition: pressure.h:69
geometry
t_geometry geometry
Definition: geometry.cpp:5
t_conv::getCounterName
const char * getCounterName(const long type) const
Definition: conv.h:344
Wind::AccelCont
realnum AccelCont
Definition: wind.h:55
colliders
ColliderList colliders(dense)
ipIRON
const int ipIRON
Definition: cddefines.h:373
ipoint
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:15
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:558
t_hcmap::lgMapDone
bool lgMapDone
Definition: hcmap.h:36
t_dense::pden
realnum pden
Definition: dense.h:108
t_yield::nelec_eject
long nelec_eject(long n, long i, long ns) const
Definition: yield.h:55
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:477
SaveGrid
void SaveGrid(FILE *pnunit, exit_type status)
Definition: save_do.cpp:4694
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
GrainVar::GrainEmission
vector< realnum > GrainEmission
Definition: grainvar.h:580
t_yield::elec_eject_frac
realnum elec_eject_frac(long n, long i, long ns, long ne) const
Definition: yield.h:48
POW2
#define POW2
Definition: cddefines.h:983
t_isoCTRL::nLyman
long int nLyman[NISO]
Definition: iso.h:352
energy
double energy(const genericState &gs)
Definition: generic_state.cpp:49
ipH1s
const int ipH1s
Definition: iso.h:29
t_hmi::HeatH2Dish_used
double HeatH2Dish_used
Definition: hmi.h:140
t_conv::nPres2Ioniz
long int nPres2Ioniz
Definition: conv.h:145
STATIC
#define STATIC
Definition: cddefines.h:118
t_save::chSaveArgs
char chSaveArgs[LIMPUN][5]
Definition: save.h:368
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:598
diatomics::lgEnabled
bool lgEnabled
Definition: h2_priv.h:352
EmissionProxy::Ploss
realnum Ploss() const
Definition: emission.h:130
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:447
ipExtraLymanLines
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
Save_Line_RT
void Save_Line_RT(FILE *ip)
Definition: save_line.cpp:328
t_save::ncSaveSkip
long int ncSaveSkip
Definition: save.h:466
t_colden::H0_21cm_lower
double H0_21cm_lower
Definition: colden.h:67
FindStrongestLineLabels
STATIC void FindStrongestLineLabels(void)
Definition: save_do.cpp:70
mole
t_mole_local mole
Definition: mole.cpp:8
t_conv::hist_pres_current
vector< double > hist_pres_current
Definition: conv.h:293
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:787
t_dense::EdenTrue
double EdenTrue
Definition: dense.h:232
pressure
t_pressure pressure
Definition: pressure.cpp:9
rfield
t_rfield rfield
Definition: rfield.cpp:9
t_ionbal::UTA_ionize_rate
double ** UTA_ionize_rate
Definition: ionbal.h:171
mean
t_mean mean
Definition: mean.cpp:16
GridGatherInCloudy
void GridGatherInCloudy(void)
Definition: grid_xspec.cpp:174
TransitionProxy::ipCont
long & ipCont() const
Definition: transition.h:489
t_conv::hist_temp_temp
vector< double > hist_temp_temp
Definition: conv.h:298
t_rfield::ConInterOut
realnum * ConInterOut
Definition: rfield.h:156
realnum
float realnum
Definition: cddefines.h:124
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:468
NUM_OUTPUT_TYPES
const int NUM_OUTPUT_TYPES
Definition: grid.h:22
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:168
t_colden::H0_21cm_upper
double H0_21cm_upper
Definition: colden.h:66
Wind::AccelTotalOutward
realnum AccelTotalOutward
Definition: wind.h:52
t_input::echo
void echo(FILE *ipOUT)
Definition: input.cpp:102
t_colden::H0_ov_Tspin
double H0_ov_Tspin
Definition: colden.h:52
cdSPEC
void cdSPEC(int Option, long int nEnergy, double ReturnedSpectrum[])
Definition: cdspec.cpp:20
EmissionProxy::Pesc
realnum & Pesc() const
Definition: emission.h:568
INPUT_LINE_LENGTH
const int INPUT_LINE_LENGTH
Definition: cddefines.h:301
t_save::nSaveEveryZone
long int nSaveEveryZone[LIMPUN]
Definition: save.h:365
max
long max(int a, long b)
Definition: cddefines.h:821
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:435
t_opac::telec
realnum telec
Definition: opacity.h:187
t_rfield::otscon
realnum * otscon
Definition: rfield.h:185
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:484
t_save::lgLineListRatio
bool lgLineListRatio[LIMPUN]
Definition: save.h:256
molecule::index
int index
Definition: mole.h:194
t_cpu_i::lgMPI_talk
bool lgMPI_talk() const
Definition: cpu.h:394
t_timesc::time_H2_Dest_here
double time_H2_Dest_here
Definition: timesc.h:48
ipRecRad
const int ipRecRad
Definition: cddefines.h:332
t_save::chLineListLabel
vector< string > chLineListLabel[LIMPUN]
Definition: save.h:252
t_hcmap::lgMapBeingDone
bool lgMapBeingDone
Definition: hcmap.h:33
t_opac::TauAbsFace
realnum * TauAbsFace
Definition: opacity.h:99
t_mesh::anu2
double anu2(size_t i) const
Definition: mesh.h:115
h2
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
GetDopplerWidth
realnum GetDopplerWidth(realnum massAMU)
Definition: temp_change.cpp:355
safe_div
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:1015
t_rfield::OccNumbIncidCont
realnum * OccNumbIncidCont
Definition: rfield.h:119
t_radius::depth_mid_zone
double depth_mid_zone
Definition: radius.h:31
iterations
t_iterations iterations
Definition: iterations.cpp:6
ipRecEsc
const int ipRecEsc
Definition: cddefines.h:328
column
double column(const genericState &gs)
Definition: generic_state.cpp:15
map_do
void map_do(FILE *io, const char *chType)
Definition: hcmap.cpp:25
t_pressure::PresGasCurr
double PresGasCurr
Definition: pressure.h:46
nWindLine
long nWindLine
Definition: cdinit.cpp:19
t_pressure::PresTotlError
double PresTotlError
Definition: pressure.h:46
optimize
t_optimize optimize
Definition: optimize.cpp:6
PrtLogLin
double PrtLogLin(double value)
Definition: save_do.cpp:724
lgMustPrintHeader
static bool lgMustPrintHeader
Definition: save_line.cpp:264
t_rfield::chLineLabel
vector< string > chLineLabel
Definition: rfield.h:212
grid
t_grid grid
Definition: grid.cpp:5
cdExecTime
double cdExecTime()
Definition: cddrive.cpp:483
radius
t_radius radius
Definition: radius.cpp:5
t_thermal::heating
double heating(long nelem, long ion)
Definition: thermal.h:186
timesc
t_timesc timesc
Definition: timesc.cpp:7
t_hmi::H2_photodissoc_used_H2g
double H2_photodissoc_used_H2g
Definition: hmi.h:119
AnuUnit
double AnuUnit(realnum energy)
Definition: service.cpp:197
t_mesh::anumin
double anumin(size_t i) const
Definition: mesh.h:139
t_save::params
SaveParams params[LIMPUN]
Definition: save.h:269
t_rfield::reflin
realnum ** reflin
Definition: rfield.h:198
t_rfield::ConEmitOut
realnum ** ConEmitOut
Definition: rfield.h:153
t_ionbal::PhotoRate_Shell
double **** PhotoRate_Shell
Definition: ionbal.h:109
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:160
t_gaunt::gauntff
double gauntff(long Z, double Te, double anu)
Definition: atmdat_gaunt.cpp:135
t_rfield::fine_opac_zone
realnum * fine_opac_zone
Definition: rfield.h:389
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:473
PrettyTransmission
realnum PrettyTransmission(long j, realnum transmission)
Definition: save_do.cpp:708
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_rfield::extin_mag_V_point
double extin_mag_V_point
Definition: rfield.h:260
chDummy
char * chDummy
Definition: save_do.cpp:815
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:76
t_hmi::UV_Cont_rel2_Draine_DB96_depth
realnum UV_Cont_rel2_Draine_DB96_depth
Definition: hmi.h:84
t_radius::Conv2PrtInten
double Conv2PrtInten
Definition: radius.h:152
ipSILICON
const int ipSILICON
Definition: cddefines.h:361
t_radius::Radius_mid_zone
double Radius_mid_zone
Definition: radius.h:31
t_rfield::ipV_filter
long int ipV_filter
Definition: rfield.h:242
ipH2p
const int ipH2p
Definition: iso.h:31
t_ionbal::CollIonRate_Ground
double *** CollIonRate_Ground
Definition: ionbal.h:118
t_save::lgSaveEveryZone
bool lgSaveEveryZone[LIMPUN]
Definition: save.h:364
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
EmissionProxy::Pdest
realnum & Pdest() const
Definition: emission.h:588
t_save::lgLinEvery
bool lgLinEvery
Definition: save.h:463
t_save::SaveHeaderDone
void SaveHeaderDone(int ipPun)
Definition: save.h:349
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:617
sprt_wl
void sprt_wl(char *chString, realnum wl)
Definition: prt.cpp:56
t_save::lgPrtIsotropicCont
bool lgPrtIsotropicCont[LIMPUN]
Definition: save.h:300
t_timesc::sound_speed_adiabatic
double sound_speed_adiabatic
Definition: timesc.h:58
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:431
t_iterations::lgLastIt
bool lgLastIt
Definition: iterations.h:47
t_save::lgFITS
bool lgFITS[LIMPUN]
Definition: save.h:377
OccupationNumberLine
double OccupationNumberLine(const TransitionProxy &t)
Definition: transition.cpp:177
t_hyperfine::Tspin21cm
double Tspin21cm
Definition: hyperfine.h:56
EmissionProxy::TauCon
realnum & TauCon() const
Definition: emission.h:498
ipH2s
const int ipH2s
Definition: iso.h:30
t_conv::lgConvTemp
bool lgConvTemp
Definition: conv.h:191
ipALUMINIUM
const int ipALUMINIUM
Definition: cddefines.h:360
mole_dominant_rates
void mole_dominant_rates(const vector< const molecule * > &debug_list, FILE *ioOut, bool lgPrintReagents, size_t NPRINT, double fprint)
Definition: save_species.cpp:589
t_phycon::EnergyBinding
double EnergyBinding
Definition: phycon.h:54
t_dynamics::Rate
double Rate
Definition: dynamics.h:77
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1069
t_save::wlLineList
vector< realnum > wlLineList[LIMPUN]
Definition: save.h:254
ChargTranPun
void ChargTranPun(FILE *ipPnunit, char *chSave)
Definition: atmdat_char_tran.cpp:1728
density
double density(const genericState &gs)
Definition: generic_state.cpp:27
t_rfield::extin_mag_V_extended
double extin_mag_V_extended
Definition: rfield.h:264
t_save::lgHashEndIter
bool lgHashEndIter[LIMPUN]
Definition: save.h:394
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:64
ExtraLymanLines
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
LIMELM
const int LIMELM
Definition: cddefines.h:307
cdColm
int cdColm(const char *chLabel, long int ion, double *theocl)
Definition: cddrive.cpp:595
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:50
t_rfield::nfine
long nfine
Definition: rfield.h:385
t_radius::drad_x_fillfac
double drad_x_fillfac
Definition: radius.h:76
Energy::Ryd
double Ryd() const
Definition: energy.h:33
molezone::den
double den
Definition: mole.h:421
t_Heavy::Valence_IP_Ryd
double Valence_IP_Ryd[LIMELM][LIMELM]
Definition: heavy.h:24
t_save::lgPrtOldStyleLogs
bool lgPrtOldStyleLogs[LIMPUN]
Definition: save.h:275
TransitionProxy::Coll
CollisionProxy Coll() const
Definition: transition.h:463
diatomics::H2_PunchLineStuff
void H2_PunchLineStuff(FILE *io, realnum xLimit, long index)
Definition: mole_h2_io.cpp:1158
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:349
generic_state.h
t_rfield::comdn
double * comdn
Definition: rfield.h:238
t_phycon::te10
double te10
Definition: phycon.h:58
t_dense::xMassDensity
realnum xMassDensity
Definition: dense.h:101
diatomics::H2_PunchDo
void H2_PunchDo(FILE *io, char chJOB[], const char chTime[], long int ipPun)
Definition: mole_h2_io.cpp:1195
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:585
dBaseStates
vector< qList > dBaseStates
Definition: taulines.cpp:16
t_hmi::H2_total
double H2_total
Definition: hmi.h:25
t_thermal::elementcool
double elementcool[LIMELM+1]
Definition: thermal.h:111
t_save::whichDiatomToPrint
diatomics * whichDiatomToPrint[LIMPUN]
Definition: save.h:307
ConvPresTempEdenIoniz
int ConvPresTempEdenIoniz(void)
Definition: conv_pres_temp_eden_ioniz.cpp:23
t_mesh::egamry
double egamry() const
Definition: mesh.h:88
t_conv::hist_temp_heat
vector< double > hist_temp_heat
Definition: conv.h:298
t_phycon::TeInit
double TeInit
Definition: phycon.h:99
t_save::ipEmisFreq
long ipEmisFreq[LIMPUN]
Definition: save.h:484
t_grid::lgOutputTypeOn
bool lgOutputTypeOn[NUM_OUTPUT_TYPES]
Definition: grid.h:66
dBaseSpecies
vector< species > dBaseSpecies
Definition: taulines.cpp:15
t_dense::eden
double eden
Definition: dense.h:201
t_opac::TauAbsGeo
realnum ** TauAbsGeo
Definition: opacity.h:90
SaveSpecies
void SaveSpecies(FILE *ioPUN, long int ipPun)
Definition: save_species.cpp:35
TransitionProxy::EnergyRyd
double EnergyRyd() const
Definition: transition.h:95
t_phycon::EnergyExcitation
double EnergyExcitation
Definition: phycon.h:47
NCHLAB
const int NCHLAB
Definition: cddefines.h:303
MAX2
#define MAX2(a, b)
Definition: cddefines.h:828
t_ionbal::RateIonizTot
double RateIonizTot(long nelem, long ion) const
Definition: ionbal.cpp:223
t_hmi::UV_Cont_rel2_Habing_TH85_depth
realnum UV_Cont_rel2_Habing_TH85_depth
Definition: hmi.h:74
lgCheckMonitors
bool lgCheckMonitors(FILE *ioMONITOR)
Definition: monitor_results.cpp:1664
EmissionProxy::damp
realnum & damp() const
Definition: emission.h:608
fprintf
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
t_phycon::te70
double te70
Definition: phycon.h:58
t_mean::xIonMean
multi_arr< double, 4 > xIonMean
Definition: mean.h:17
t_Heavy::chShell
char chShell[7][3]
Definition: heavy.h:31
t_save::chSpeciesDominantRates
string chSpeciesDominantRates[LIMPUN]
Definition: save.h:480
t_save::chConSavEnr
const char * chConSavEnr[LIMPUN]
Definition: save.h:387
t_mole_global::list
MoleculeList list
Definition: mole.h:365
saveFITSfile
void saveFITSfile(FILE *io, int option)
Definition: save_fits.cpp:83
t_DoppVel::lgTurb_pressure
bool lgTurb_pressure
Definition: doppvel.h:42
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:33
prt_line_inlist
void prt_line_inlist(FILE *ioOUT, const char *label, realnum wvlng)
Definition: prt.cpp:173
t_dynamics::Source
double ** Source
Definition: dynamics.h:80
AGN_Hemis
STATIC void AGN_Hemis(FILE *ioPUN)
Definition: save_do.cpp:4467
SaveSpeciesPseudoCont
void SaveSpeciesPseudoCont(const long ipPun, const string &speciesLabel)
Definition: species_pseudo_cont.cpp:403
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:1006
t_save::chSaveSpecies
vector< string > chSaveSpecies[LIMPUN]
Definition: save.h:370
t_dynamics::time_elapsed
double time_elapsed
Definition: dynamics.h:105
t_phycon::telogn
double telogn[7]
Definition: phycon.h:86
t_conv::ntypes
size_t ntypes(void) const
Definition: conv.h:307
ipCARBON
const int ipCARBON
Definition: cddefines.h:353
t_rfield::TotDiff2Pht
realnum * TotDiff2Pht
Definition: rfield.h:179
t_grid::seqNum
long seqNum
Definition: grid.h:72
DynaPunchTimeDep
void DynaPunchTimeDep(FILE *ipPnunit, const char *chJob)
Definition: dynamics.cpp:2057
save_average
void save_average(long int ipPun)
Definition: save_average.cpp:128
t_iso_sp::numLevels_max
long int numLevels_max
Definition: iso.h:524
hcmap
t_hcmap hcmap
Definition: hcmap.cpp:23
GammaPrt
void GammaPrt(long int ipLoEnr, long int ipHiEnr, long int ipOpac, FILE *ioFILE, double total, double threshold)
Definition: cont_gammas.cpp:244
SaveLineIntensity
STATIC void SaveLineIntensity(FILE *ioPUN, long int ipPun, realnum Threshold)
Definition: save_do.cpp:4150
t_LineSave::nsum
long int nsum
Definition: lines.h:87
t_pressure::pinzon_PresIntegElecThin
realnum pinzon_PresIntegElecThin
Definition: pressure.h:75
t_struc::testr
realnum * testr
Definition: struc.h:25
GrainVar::SilicateEmission
vector< realnum > SilicateEmission
Definition: grainvar.h:582
t_mesh::anumax
double anumax(size_t i) const
Definition: mesh.h:143
t_pressure::PresTurbCurr
double PresTurbCurr
Definition: pressure.h:42
fixit
#define fixit(a)
Definition: cddefines.h:416
t_rfield::ConEmitReflec
realnum ** ConEmitReflec
Definition: rfield.h:147
gv
GrainVar gv
Definition: grainvar.cpp:5
hmi
t_hmi hmi
Definition: hmi.cpp:5
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
cpu
static t_cpu cpu
Definition: cpu.h:423
POW3
#define POW3
Definition: cddefines.h:990
ShowMe
void ShowMe(void)
Definition: service.cpp:205
dBaseTrans
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
save
t_save save
Definition: save.cpp:5
t_phycon::te
double te
Definition: phycon.h:21
t_rfield::chCumuType
char chCumuType[5]
Definition: rfield.h:333
t_thermal::dCooldT
double dCooldT
Definition: thermal.h:139
plankf
double plankf(long int ip)
Definition: service.cpp:1869
t_save::Resolution
realnum Resolution
Definition: save.h:475
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:348
SaveResults
STATIC void SaveResults(FILE *ioPUN)
Definition: save_do.cpp:4536
SaveNewContinuum
STATIC void SaveNewContinuum(FILE *ioPUN)
Definition: save_do.cpp:4374
t_rfield::nflux
long int nflux
Definition: rfield.h:48
t_save::emisfreq
Energy emisfreq[LIMPUN]
Definition: save.h:483
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:668
flux_correct_isotropic
double flux_correct_isotropic(const bool lgSaveIsotr, const int nEmType, const int iflux)
Definition: rfield.cpp:112
t_colden::colden
realnum colden[NCOLD]
Definition: colden.h:32
t_save::lgSaveHeader
bool lgSaveHeader(int ipPun) const
Definition: save.h:345
t_rfield::ConRefIncid
realnum ** ConRefIncid
Definition: rfield.h:159
t_save::ipConPun
long int ipConPun
Definition: save.h:390
LinSv
Definition: lines.h:157
t_save::nLineList
long nLineList[LIMPUN]
Definition: save.h:250
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:529
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:519
t_save::chSave
char chSave[LIMPUN][5]
Definition: save.h:306
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:468
t_mesh::getResolutionScaleFactor
double getResolutionScaleFactor() const
Definition: mesh.h:92
Wind::windv
realnum windv
Definition: wind.h:18
ipLine
static long int * ipLine
Definition: prt_linesum.cpp:14
CollisionProxy::ColUL
double ColUL(const ColliderList &colls) const
Definition: collision.h:106
EmissionProxy::TauIn
realnum & TauIn() const
Definition: emission.h:458
t_warnings::lgWarngs
bool lgWarngs
Definition: warnings.h:44
t_LineSave::ipass
long int ipass
Definition: lines.h:96
called
t_called called
Definition: called.cpp:4
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:363
t_ionbal::CotaRate
realnum CotaRate[LIMELM]
Definition: ionbal.h:236
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
GrainVar::rate_h2_form_grains_used_total
double rate_h2_form_grains_used_total
Definition: grainvar.h:576
t_cpu_i::chExitStatus
const string & chExitStatus(exit_type s) const
Definition: cpu.h:399
ipSODIUM
const int ipSODIUM
Definition: cddefines.h:358
t_ionbal::RR_rate_coef_used
double ** RR_rate_coef_used
Definition: ionbal.h:207
lgPopsFirstCall
static bool lgPopsFirstCall
Definition: save_do.cpp:4198
t_thermal::ctot
double ctot
Definition: thermal.h:130
Wind::fmul
realnum fmul
Definition: wind.h:65
Save1Line
void Save1Line(const TransitionProxy &t, FILE *io, realnum xLimit, long index, realnum DopplerWidth)
Definition: save_do.cpp:4316
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:166
EmissionProxy::pump
double & pump() const
Definition: emission.h:518
Generated by   doxygen 1.8.5