cloudy  trunk
save_do.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*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 /*SaveFeII_cont return appropriate FeII banded continuum, inward, outward, or total */
18 /*FindStrongestLineLabels find strongest lines contributing to point in continuum energy mesh, output in some save commands */
19 #include "cddefines.h"
20 #include "cddrive.h"
21 #include "physconst.h"
22 #include "mean.h"
23 #include "taulines.h"
24 #include "struc.h"
25 #include "iso.h"
26 #include "mole.h"
27 #include "hyperfine.h"
28 #include "rt.h"
29 #include "lines_service.h"
30 #include "doppvel.h"
31 #include "dense.h"
32 #include "h2.h"
33 #include "magnetic.h"
34 #include "hydrogenic.h"
35 #include "secondaries.h"
36 #include "grainvar.h"
37 #include "lines.h"
38 #include "dynamics.h"
39 #include "colden.h"
40 #include "continuum.h"
41 #include "ionbal.h"
42 #include "yield.h"
43 #include "prt.h"
44 #include "iterations.h"
45 #include "heavy.h"
46 #include "conv.h"
47 #include "geometry.h"
48 #include "called.h"
49 #include "helike.h"
50 #include "opacity.h"
51 #include "rfield.h"
52 #include "phycon.h"
53 #include "timesc.h"
54 #include "radius.h"
55 #include "atomfeii.h"
56 #include "monitor_results.h"
57 #include "thermal.h"
58 #include "wind.h"
59 #include "hmi.h"
60 #include "pressure.h"
61 #include "elementnames.h"
62 #include "ipoint.h"
63 #include "gammas.h"
64 #include "atmdat.h"
65 #include "hcmap.h"
66 #include "input.h"
67 #include "save.h"
68 #include "optimize.h"
69 #include "warnings.h"
70 #include "grid.h"
71 #include "mole_priv.h"
72 #include "atomfeii.h"
73 
74 // this needs C linkage since it is passed as a pointer to a C library routine
75 extern "C" int wavelength_compare (const void * a, const void * b)
76 {
77  LinSv a1 = *(LinSv*)a;
78  LinSv b1 = *(LinSv*)b;
79  /* comparison is b-a so we get inverse wavelength order (increasing energy order) */
80  if( b1.wavelength > a1.wavelength )
81  return 1;
82  else if( b1.wavelength < a1.wavelength )
83  return -1;
84  else
85  return 0;
86 }
87 
88 // find strongest lines contributing to point in continuum energy mesh, output in some save commands
89 STATIC void FindStrongestLineLabels( void )
90 {
91  long low_index=0;
92  long high_index=0;
93  long j_min = 0;
94  double MaxFlux = 0.;
95  long ipMaxFlux = 0;
96  long j = 0;
97 
98  ASSERT( LineSave.ipass==1 );
99  // copy and sort
100  memcpy(LineSvSortWL, LineSv, (unsigned)LineSave.nsum*sizeof( LinSv ));
101  qsort((void *)LineSvSortWL,(size_t)LineSave.nsum,sizeof( LinSv ),wavelength_compare);
102 
103  while( rfield.anu[j_min]+0.5*rfield.widflx[j_min] < RYDLAM/LineSvSortWL[0].wavelength)
104  j_min++;
105 
106  for( j=0; j<rfield.nflux; j++ )
107  {
108  if( j < j_min )
109  {
110  strcpy( rfield.chLineLabel[j], " " );
111  continue;
112  }
113 
114  ASSERT( LineSvSortWL[low_index].wavelength != 0. );
115 
116  while( RYDLAM/LineSvSortWL[low_index].wavelength < rfield.anu[j]-0.5*rfield.widflx[j] && low_index < LineSave.nsum-1 )
117  {
118  low_index++;
119  if( LineSvSortWL[low_index].wavelength == 0. )
120  {
121  // hit the end of real wavelengths. Pad rest of labels with spaces
122  for( long j1=j; j1<rfield.nflux; j1++ )
123  strcpy( rfield.chLineLabel[j1], " " );
124  return;
125  }
126  }
127 
128  high_index = low_index;
129  ASSERT( LineSvSortWL[high_index].wavelength != 0. );
130 
131  while( RYDLAM/LineSvSortWL[high_index].wavelength < rfield.anu[j]+0.5*rfield.widflx[j] && high_index < LineSave.nsum-1 )
132  {
133  high_index++;
134  if( LineSvSortWL[high_index].wavelength == 0. )
135  {
136  high_index--;
137  break;
138  }
139  }
140  // while loop found first one greater than j bin, decrement again to get back into j bin
141  high_index--;
142 
143  ASSERT( LineSvSortWL[low_index].wavelength > 0. );
144  ASSERT( LineSvSortWL[high_index].wavelength > 0. );
145  ASSERT( RYDLAM/LineSvSortWL[low_index].wavelength >= rfield.anu[j]-0.5*rfield.widflx[j] );
146  ASSERT( RYDLAM/LineSvSortWL[high_index].wavelength <= rfield.anu[j]+0.5*rfield.widflx[j] );
147 
148  MaxFlux = 0.;
149  ipMaxFlux = 0;
150 
151  for( long k = low_index; k <= high_index; k++ )
152  {
153  if( strcmp( LineSvSortWL[k].chALab, "Coll" )==0 ||
154  strcmp( LineSvSortWL[k].chALab, "Heat" )==0 ||
155  strcmp( LineSvSortWL[k].chALab, "Pump" )==0 ||
156  strcmp( LineSvSortWL[k].chALab, "pump" )==0 ||
157  strcmp( LineSvSortWL[k].chALab, "nInu" )==0 ||
158  strcmp( LineSvSortWL[k].chALab, "nFnu" )==0 ||
159  strcmp( LineSvSortWL[k].chALab, "InwT" )==0 ||
160  strcmp( LineSvSortWL[k].chALab, "InwC" )==0 ||
161  strcmp( LineSvSortWL[k].chALab, "Inwd" )==0 ||
162  strcmp( LineSvSortWL[k].chALab, "Ca A" )==0 ||
163  strcmp( LineSvSortWL[k].chALab, "Ca B" )==0 ||
164  strcmp( LineSvSortWL[k].chALab, "Pho+" )==0 ||
165  strcmp( LineSvSortWL[k].chALab, "Pcon" )==0 ||
166  strcmp( LineSvSortWL[k].chALab, "Q(H)" )==0 ||
167  strcmp( LineSvSortWL[k].chALab, "Unit" )==0 )
168  continue;
169 
170  if( LineSvSortWL[k].SumLine[0] > MaxFlux )
171  {
172  MaxFlux = LineSvSortWL[k].SumLine[0];
173  ipMaxFlux = k;
174  }
175  }
176 
177  /* line label */
178  if( ipMaxFlux > 0 )
179  strcpy( rfield.chLineLabel[j], LineSvSortWL[ipMaxFlux].chALab );
180  }
181 
182  return;
183 }
184 
185 # ifdef USE_NLTE7
186 //Run "Save NLTE"
187 static void runNLTE(int ipPun)
188 {
189  long nelem = ipNEON;
190  int nouter[11] = {0,2,3,3,3,4,3,5,5,6,0};
191  ASSERT(dense.gas_phase[nelem] != 0.);
192  // Section 1
193  fprintf( save.ipPnunit[ipPun], "data Ferland University of Kentucky Cloudy 10 11-5-2011\n");
194  fprintf( save.ipPnunit[ipPun], "case NeXY6\n");
195  fprintf( save.ipPnunit[ipPun], "code Cloudy\n");
196  fprintf( save.ipPnunit[ipPun], "atom Ne 10\n");
197  fprintf( save.ipPnunit[ipPun], "calctime %11.4e %11.4e\n",4.5,30.);
198  //Section 2
200  //Avg Charge
201  double avgchrg = 0.;
202  double m2 = 0.;
203  double m3 = 0.;
204  //double m4 = 0.;
205  for( int ion = 1; ion < nelem-1; ion++)
206  {
207  avgchrg += (dense.xIonDense[nelem][ion]/dense.gas_phase[nelem] * ion);
208  }
209  //Central Moments of Charge Distribution
210  for( int ion = 1; ion < nelem-1; ion++)
211  {
212  m2 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW2((ion - avgchrg));
213  m3 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW3((ion - avgchrg));
214  //m4 += dense.xIonDense[nelem][ion]/dense.gas_phase[nelem]*POW4((ion - avgchrg));
215  }
216  //Specific internal energy
217  double Eint = 0.;
218  double partfun = 0.;
219  long count_nrglvls = 0;
220  for( int ion = 1; ion < nelem-1; ion++)
221  {
222  // Find ipSpecies
223  long ipSpec = 0;
224  for( int i = 0; i < nSpecies; i++)
225  {
226  ipSpec = i;
227  if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion)
228  break;
229  ipSpec = -1;
230  }
231 
232  if( ipSpec != -1)
233  {//This requires a loop over energy levels
234  for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++)
235  {
236 
237  Eint += dBaseStates[ipSpec][lvl].Pop()*dBaseStates[ipSpec][lvl].energy().eV();
238  partfun += dBaseStates[ipSpec][lvl].g()*exp(-1*dBaseStates[ipSpec][lvl].energy().eV()/
239  phycon.te_eV);
240  count_nrglvls += 1;
241  }
242  }
243  else
244  {
245  printf("\nNot using :\t%li\t%i\n",nelem+1,ion+1);
246  //cdEXIT(EXIT_FAILURE);
247  }
248  }
249  Eint = Eint / dense.gas_phase[nelem];
250 
251 
252  fprintf( save.ipPnunit[ipPun],"\nsummary_quantities\n");
253  fprintf( save.ipPnunit[ipPun],"plasma %11.4e %11.4e\n",phycon.te_eV,dense.eden);
254  fprintf( save.ipPnunit[ipPun],"time %11.4e\n",0.);
255  fprintf( save.ipPnunit[ipPun],"zbar %11.4e\n",avgchrg);
256  fprintf( save.ipPnunit[ipPun],"m2 %11.4e\n",m2);
257  fprintf( save.ipPnunit[ipPun],"m3 %11.4e\n",m3);
258  //fprintf( save.ipPnunit[ipPun],"m4 %11.3e\n",m4);
259  fprintf( save.ipPnunit[ipPun],"eint %11.4e\n",0.); //%11.3e\n",Eint);
260  fprintf( save.ipPnunit[ipPun],"deintdt %11.4e\n",0.); //%11.4e\n",thermal.dCooldT);
261  fprintf( save.ipPnunit[ipPun],"pfn %11.4e\n",partfun);
262  fprintf( save.ipPnunit[ipPun],"nmax_eff %i\n",0);
263  fprintf( save.ipPnunit[ipPun],"ploss %11.4e %11.4e %11.4e %9.3e\n",0.,0.,0.,thermal.totcol);
264  //Section 3
265  fprintf( save.ipPnunit[ipPun],"\nion_stages %li\n",nelem-2);
266 
267  //Photoionization Rate
268  double PIR[nelem+1];
269  //autoionization Rate
270  double autoion[nelem+1];
271  //This requires a loop over ion stages
272  for( int ion = 1; ion < nelem-1; ion++)
273  {
274  double tempRecomb = 0.;
275  double f_aauto = 0.;
276  double f_aphoto = 0.;
277  double f_acoll = 0.;
278  double f_Scoll = 0.;
279  double f_Sphoto = 0.;
280  double f_auto = 0.;
281  autoion[ion] = 0.;
282  PIR[ion] = 0.;
283  //This deals with the fact that there are no recombinations FROM atoms
284  if( ion != 0 && ionbal.RateRecomTot[nelem][ion-1] != 0.)
285  {
286  tempRecomb = ionbal.RateRecomTot[nelem][ion-1];
287  f_aauto = dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion-1]/tempRecomb;
288  f_aphoto = dense.eden*ionbal.RR_rate_coef_used[nelem][ion-1]/tempRecomb;
289  f_acoll = dense.eden*ionbal.CotaRate[ion-1]/tempRecomb;
290  }
291  //Photoionization rate
292  for( int i=0; i < Heavy.nsShells[nelem][ion]; i++)
293  {
294  //Add up the photoionization rates for all of the shells
295  if( ion != nelem+1)
296  {
297  PIR[ion] += ionbal.PhotoRate_Shell[nelem][ion][i][0];
298  }
299  }
300  //Compute the fractional ionizations
301  if( ion != nelem+1 && ionbal.RateIonizTot(nelem,ion) != 0.)
302  {
303  f_Scoll = ionbal.CollIonRate_Ground[nelem][ion][0]/ionbal.RateIonizTot(nelem,ion);
304  f_Sphoto = PIR[ion]/ionbal.RateIonizTot(nelem,ion);
305  f_auto = autoion[ion]/ionbal.RateIonizTot(nelem,ion);
306  }
307  //Autoionizations
308  autoion[ion] = ionbal.UTA_ionize_rate[nelem][ion] + secondaries.csupra[nelem][ion];
309  if( ion > 0)
310  {
311  autoion[ion] += 0.0; //auger[ion-1];
312  }
313  fprintf( save.ipPnunit[ipPun],"ion %li %11.4e %i\n"
314  " %11.3e %11.3e %11.3e %11.3e\n"
315  " %11.3e %11.3e %11.3e %11.3e\n\n",
316  nelem+1-ion,
317  dense.xIonDense[nelem][ion]/dense.gas_phase[nelem],
318  nouter[ion],
319  ionbal.RateIonizTot(nelem,ion),
320  f_Scoll,
321  f_Sphoto,
322  f_auto,
323  tempRecomb,
324  f_acoll,
325  f_aphoto,
326  f_aauto);
327  }
328 
329  //Section 4
330  fprintf( save.ipPnunit[ipPun],"\nenergy_levels %li\n",count_nrglvls);
331  //Collisional Destruction Rate Bound-Free - Collisional Ionization
332  double GcollBF = 0.0;
333  //Photo Dest Rate Bound-Free - Photoionization
334  double GphotoBF = 0.0;
335  //Autoionization Dest Rate
336  double GautoBF = 0.0;
337  // Total Dest Rate
338  double Gtotal = 0.0;
339  //Collisional Creation Rate Bound-Free - Collisional Ionization
340  double QcollBF = 0.0;
341  //Photo Creation Rate Bound-Free - Photoionization
342  double QphotoBF = 0.0;
343  //Autoionization Creation Rate
344  double QautoBF = 0.0;
345  // Total Creation Rate
346  double Qtotal = 0.0;
347  //total pop of all levels
348  double totallevelpop = 0.;
349 
350  for( int ion = 1; ion < nelem-1; ion++)
351  {
352  // Find ipSpecies
353  long ipSpec = 0;
354  for( int i = 0; i < nSpecies; i++)
355  {
356  ipSpec = i;
357  if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion)
358  break;
359  ipSpec = -1;
360  }
361  if( ipSpec != -1)
362  {
363  for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++)
364  {
365  //total level population
366  totallevelpop += dBaseStates[ipSpec][lvl].Pop();
367  }
368  }
369  }
370 
371  for( int ion = 1; ion < nelem-1; ion++)
372  {
373  // Find ipSpecies
374  long ipSpec = 0;
375  for( int i = 0; i < nSpecies; i++)
376  {
377  ipSpec = i;
378  if(dBaseStates[i][0].nelem()-1 == nelem && dBaseStates[i][0].IonStg()-1 == ion)
379  break;
380  ipSpec = -1;
381  }
382  //Energy Level Correction Factor - Scale energy levels to ground state of atomic Neon
383  //using ionization potentials
384  double ELCF = 0.;
385  for( int i = 0; i < ion; i++ )
386  {
387  ELCF += Heavy.Valence_IP_Ryd[nelem][i]*EVRYD;
388  }
389 
390  if( ipSpec != -1)
391  {
392 
393  //This requires a loop over energy levels
394  for( int lvl = 0; lvl < dBaseSpecies[ipSpec].numLevels_max; lvl++)
395  {
396 
397  //Give bound-free Destruction rates for the ground state, zero otherwise
398  if( lvl == 0 && ion != nelem+1)
399  {
400  GcollBF = ionbal.CollIonRate_Ground[nelem][ion][0];
401  GphotoBF = PIR[ion];
402  GautoBF = autoion[ion];
403  }
404  else
405  {
406  GcollBF = 0.;
407  GphotoBF = 0.;
408  GautoBF = 0.;
409  }
410  //Find the Destruction total rate
411  Gtotal = dBaseStates[ipSpec][lvl].DestCollBB() +
412  dBaseStates[ipSpec][lvl].DestPhotoBB() +
413  GcollBF + GphotoBF + GautoBF;
414  //Store Dest Rate Fractions
415  double f_GcollBB,
416  f_GphotoBB,
417  f_GcollBF,
418  f_GphotoBF,
419  f_GautoBF;
420 
421  if( Gtotal != 0.)
422  {
423  f_GcollBB = dBaseStates[ipSpec][lvl].DestCollBB()/Gtotal;
424  f_GphotoBB = dBaseStates[ipSpec][lvl].DestPhotoBB()/Gtotal;
425  f_GcollBF = GcollBF/Gtotal;
426  f_GphotoBF = GphotoBF/Gtotal;
427  f_GautoBF = GautoBF/Gtotal;
428  }
429  else
430  {
431  f_GcollBB = 0.;
432  f_GphotoBB = 0.;
433  f_GcollBF = 0.;
434  f_GphotoBF = 0.;
435  f_GautoBF = 0.;
436  }
437  //Give bound-free Creation rates for the ground state, zero otherwise
438  if( lvl == 0 && ion != 0 )
439  {
440  QautoBF = dense.eden*ionbal.DR_Badnell_rate_coef[nelem][ion-1];
441  QphotoBF = dense.eden*ionbal.RR_rate_coef_used[nelem][ion-1];
442  QcollBF = dense.eden*ionbal.CotaRate[ion-1];
443  }
444  else
445  {
446  QautoBF = 0.;
447  QphotoBF = 0.;
448  QcollBF = 0.;
449  }
450 
451  //Find the Creation total rate
452  Qtotal = dBaseStates[ipSpec][lvl].CreatCollBB() +
453  0.0 + QcollBF + QphotoBF + QautoBF;
454  //Get Creation Rate Fractions
455  double f_QcollBB,
456  f_QphotoBB,
457  f_QcollBF,
458  f_QphotoBF,
459  f_QautoBF;
460 
461  if( Qtotal != 0.)
462  {
463  f_QcollBB = dBaseStates[ipSpec][lvl].CreatCollBB()/Qtotal;
464  f_QphotoBB = 0.0/Qtotal;
465  f_QcollBF = QcollBF/Qtotal;
466  f_QphotoBF = QphotoBF/Qtotal;
467  f_QautoBF = QautoBF/Qtotal;
468  }
469  else
470  {
471  f_QcollBB = 0.;
472  f_QphotoBB = 0.;
473  f_QcollBF = 0.;
474  f_QphotoBF = 0.;
475  f_QautoBF = 0.;
476  }
477 
478  //Section 4 print
479  fprintf( save.ipPnunit[ipPun],"elev %li %3i %11.4e %11.6e %11.4e\n"
480  " %11.3e %11.3e %11.3e %11.3e %11.3e %11.3e\n"
481  " %11.3e %11.3e %11.3e %11.3e %11.3e %11.3e\n"
482  " %i %i %i\n\n",
483  nelem+1-ion,
484  lvl+1,
485  dBaseStates[ipSpec][lvl].g(),
486  dBaseStates[ipSpec][lvl].energy().eV()+ELCF,
487  dBaseStates[ipSpec][lvl].Pop()/totallevelpop,
488  Gtotal*dBaseStates[ipSpec][lvl].Pop(),
489  f_GcollBB,
490  f_GphotoBB,
491  f_GcollBF,
492  f_GphotoBF,
493  f_GautoBF,
494  Qtotal,
495  f_QcollBB,
496  f_QphotoBB,
497  f_QcollBF,
498  f_QphotoBF,
499  f_QautoBF,
500  0,
501  0,
502  nouter[ion]);
503  }
504  }
505  else
506  printf("\nNot using :\t%li\t%i for levels\n",nelem+1,ion+1);
507  }
508  return;
509 }
510 # endif
511 
512 
513 // implements the absorption option on the
514 // set save line width command
515 inline realnum PrettyTransmission(long j, realnum transmission)
516 {
517  if( save.ResolutionAbs < realnum(0.) )
518  // option to conserve energy
519  return transmission;
520  else
521  {
522  realnum corr = save.ResolutionAbs*rfield.widflx[j]/rfield.anu[j];
523  return realnum(max(0., 1. - (1.-transmission)*corr ));
524  }
525 }
526 
527 /*SaveResults1Line do single line of output for the save results and save line intensity commands */
528 /* the number of emission lines across one line of printout */
529 STATIC void SaveResults1Line(
530  FILE * ioPUN,
531  /* 4 char + null string */
532  const char *chLab,
533  realnum wl,
534  double xInten,
535  const char *chFunction);/* null terminated string saying what to do */
536 
537 /*SaveGaunts called by save gaunts command to output gaunt factors */
538 STATIC void SaveGaunts(FILE* ioPUN);
539 
540 /*SaveResults save results from save results command */
541 /*SaveResults1Line do single line of output for the save results and save line intensity commands */
542 STATIC void SaveResults(FILE* ioPUN);
543 
544 STATIC void SaveLineStuff(
545  FILE * ioPUN,
546  const char *chJob ,
547  realnum xLimit);
548 
549 /* save h emission, for chapter 4, routine is below */
550 STATIC void AGN_Hemis(FILE *ioPUN );
551 
552 /*SaveNewContinuum produce the 'save new continuum' output */
553 STATIC void SaveNewContinuum(FILE * ioPUN );
554 
555 /*SaveLineIntensity produce the 'save lines intensity' output */
556 STATIC void SaveLineIntensity(FILE * ioPUN , long int ipPun, realnum Threshold);
557 
558 char *chDummy;
559 
560 /*SaveFeII_cont return appropriate FeII banded continuum, inward, outward, or total */
561 STATIC realnum SaveFeII_cont( long int ipCont , long ipFeII_Cont_type )
562 {
563 
564  DEBUG_ENTRY( "SaveFeII_cont()" );
565 
566  // [0]=wavelength, [1] inward; [2] outward, 3=> total
567  if( ipFeII_Cont_type == 3 )
568  return( FeII_Cont[ipCont][1]+FeII_Cont[ipCont][2] );
569  else
570  return(FeII_Cont[ipCont][ipFeII_Cont_type]);
571 }
572 
573 void SaveDo(
574  /* chTime is null terminated 4 char string, either "MIDL" or "LAST" */
575  const char *chTime)
576 {
577  bool lgTlkSav;
578  long int
579  ipPun,
580  i,
581  i1,
582  ion,
583  ipConMax,
584  ipHi,
585  ipLinMax,
586  ipLo,
587  ips,
588  j,
589  jj,
590  limit,
591  nelem;
592  double ConMax,
593  RateInter,
594  av,
595  conem,
596  eps,
597  flxatt,
598  flxcor,
599  flxin,
600  flxref,
601  flxtrn,
602  fout,
603  fref,
604  fsum,
605  opConSum,
606  opLinSum,
607  stage,
608  sum,
609  texc,
610  xLinMax;
611 
612  DEBUG_ENTRY( "SaveDo()" );
613 
614  /*
615  * the "last" option on save command, to save on last iteration,
616  * is parsed at the top of the loop in only one place.
617  * no further action is needed at all for save last to work
618  * ok throughout this routine
619  */
620 
621  /*
622  * each branch can have a test whether chTime is or is not "LAST"
623  *
624  * if( strcmp(chTime,"LAST") == 0 ) <== print after iteration is complete
625  *
626  * if "LAST" then this is last call to routine after iteration complete
627  * save only if "LAST" when results at end of iteration are needed
628  *
629  * if( strcmp(chTime,"LAST") != 0 ) <== print for every zone
630  *
631  * test for .not."LAST" is for every zone result, where you do not
632  * want to save last zone twice
633  */
634 
635  /* return if no save to do */
636  if( save.nsave < 1 )
637  {
638  return;
639  }
640 
641  /* during a grid calculation this routine saves grid points after
642  * cloudy is called. we may output it below */
643  if( grid.lgGrid )
644  {
645  if( chTime[0]=='L' )
646  GridGatherInCloudy();
647  }
648 
649  // sort line labels if this is last call, this avoids multiple calls if several
650  // output options need sorted labels and is safer since labels will be sorted in
651  // case new code is added that reports the strong lines. The disadvantage is that
652  // we sort even if the labels are not used
653  if( strcmp(chTime,"LAST") == 0 )
654  {
655  // sort emission line intensities so strongest lines are reported
656  FindStrongestLineLabels();
657  }
658 
659  for( ipPun=0; ipPun < save.nsave; ipPun++ )
660  {
661  if( save.lgPunHeader[ipPun] && !nMatch(save.chHeader[ipPun],save.chNONSENSE) )
662  {
663  fprintf( save.ipPnunit[ipPun], "%s", save.chHeader[ipPun] );
664  save.lgPunHeader[ipPun] = false;
665  }
666 
667  /* this global variable to remember where in the save stack we are */
668  save.ipConPun = ipPun;
669 
670  /* used to identify case where no key found */
671  bool lgNoHitFirstBranch = false;
672 
673  /* iterations.lgLastIt is true if this is last iteration
674  * lgPunLstIter set true if 'last' key occurred on save command
675  * normally is false. This will skip saving if last set and
676  * this is not last iteration */
677  /* IMPORTANT: there is a second, identical if-statement halfway
678  * down this routine. Any changes here should be copied there! */
679  if( iterations.lgLastIt || !save.lgPunLstIter[ipPun] ||
680  // if the sim aborted, make sure that the punch output is still done.
681  // for MIDL output this is not very useful as all previous zones from
682  // the iteration where the abort occured will be missing, but for LAST
683  // output it is important to print at least placeholders in grid runs.
684  ( lgAbort && strcmp(chTime,"LAST") == 0 ) ||
685  ( dynamics.lgTimeDependentStatic && dynamics.lgStatic_completed ) )
686  {
687 
688  if( strcmp(save.chSave[ipPun],"ABUN") == 0 )
689  {
690  /* save abundances vs depth */
691  if( strcmp(chTime,"LAST") != 0 )
692  {
693  fprintf( save.ipPnunit[ipPun], "%.2f",
694  log10(MAX2(SMALLFLOAT,dense.gas_phase[ipHYDROGEN])) );
695  for( nelem=ipHELIUM; nelem < LIMELM; nelem++ )
696  {
697  /* >>chng 05 feb 03, protect against non-positive abundances,
698  * bug caught by Marcelo Castellanos */
699  fprintf( save.ipPnunit[ipPun], "\t%.2f",
700  log10(MAX2(SMALLFLOAT,dense.gas_phase[nelem])) );
701  }
702  fprintf( save.ipPnunit[ipPun], "\n" );
703  }
704  }
705 
706  else if( strcmp(save.chSave[ipPun],"21CM") == 0 )
707  {
708  /* save information about 21 cm line */
709  if( strcmp(chTime,"LAST") != 0 )
710  {
711  fprintf( save.ipPnunit[ipPun],
712  "%.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",
713  /* depth, cm */
714  radius.depth_mid_zone,
715  hyperfine.Tspin21cm ,
716  phycon.te ,
717  /* temperature from Lya - 21 cm optical depth ratio */
718  3.84e-7* iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon() /
719  SDIV( HFLines[0].Emis().TauCon() ),
720  /*TexcLine( &iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ),*/
721  (*HFLines[0].Lo()).Pop() ,
722  (*HFLines[0].Hi()).Pop() ,
723  OccupationNumberLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ),
724  HFLines[0].Emis().TauCon() ,
725  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauCon(),
726  HFLines[0].Emis().PopOpc(),
727  /* term in () is density (cm-3) of 1s, so this is n(1s) / Ts */
728  (iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop())/SDIV( hyperfine.Tspin21cm),
729  /* why was above multiplied by this following term? */
730  /* *HFLines[0].EnergyErg/BOLTZMANN/4.,*/
731  HFLines[0].Emis().TauIn(),
732  TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) ) ,
733  colden.H0_ov_Tspin,
734  /*>>chng 27 mar, GS, integrated 21cm spin temperature*/
735  colden.H0_21cm_lower,
736  colden.H0_21cm_upper,
737  -0.068/log((colden.H0_21cm_upper/3.)/colden.H0_21cm_lower)
738  );
739  }
740  }
741 
742  else if( strcmp(save.chSave[ipPun],"AGES") == 0 )
743  {
744  /* save timescales vs depth */
745  if( strcmp(chTime,"LAST") != 0 )
746  {
747  int ipCO, ipOH;
748  ipCO = findspecies("CO")->index;
749  ipOH = findspecies("OH")->index;
750  fprintf( save.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
751  /* depth, cm */
752  radius.depth_mid_zone,
753  /* cooling timescale */
754  dense.pden*BOLTZMANN*1.5*phycon.te/ thermal.htot,
755  /* H2 destruction timescale */
756  timesc.time_H2_Dest_here,
757  /* CO destruction timescale */
758  1./SDIV((ipCO != -1) ? mole.species[ipCO].snk : 0.),
759  /* OH destruction timescale */
760  1./SDIV((ipOH != -1) ? mole.species[ipOH].snk : 0.),
761  /* H recombination timescale */
762  1./(dense.eden*2.90e-10/(phycon.te70*phycon.te10/phycon.te03)) );
763  }
764  }
765 
766  else if( strcmp(save.chSave[ipPun]," AGN") == 0 )
767  {
768  if( strcmp(chTime,"LAST") == 0 )
769  {
770  if( strcmp( save.chSaveArgs[ipPun], "HECS" ) == 0 )
771  {
772  /* this routine is in helike.c */
773  AGN_He1_CS(save.ipPnunit[ipPun]);
774  }
775  if( strcmp( save.chSaveArgs[ipPun], "HEMI" ) == 0 )
776  {
777  /* save h emiss, for chapter 4, routine is below */
778  AGN_Hemis(save.ipPnunit[ipPun]);
779  }
780  else
781  {
782  fprintf( ioQQQ, " SaveDo does not recognize flag %4.4s set for AGN save. This is impossible.\n",
783  save.chSave[ipPun] );
784  ShowMe();
785  cdEXIT(EXIT_FAILURE);
786  }
787  }
788  }
789 
790  else if( strcmp(save.chSave[ipPun],"MONI") == 0 )
791  {
792  if( strcmp(chTime,"LAST") == 0 )
793  {
794  /* save the monitor output */
795  lgCheckMonitors(save.ipPnunit[ipPun]);
796  }
797  }
798 
799  else if( strcmp(save.chSave[ipPun],"AVER") == 0 )
800  {
801  if( strcmp(chTime,"LAST") == 0 )
802  {
803  /* save the averages output */
804  save_average( ipPun );
805  }
806  }
807 
808  else if( strncmp(save.chSave[ipPun],"CHA",3) == 0 )
809  {
810  if( strcmp(chTime,"LAST") == 0 )
811  {
812  /* one of the charge transfer options, all in chargtran.c */
813  ChargTranPun( save.ipPnunit[ipPun] , save.chSave[ipPun] );
814  }
815  }
816 
817  else if( strcmp( save.chSave[ipPun],"CHIA") == 0)
818  {
819  static bool lgRunOnce = true;
820  if( lgRunOnce )
821  {
822  lgRunOnce = false;
823  // save chianti collision data in physical units
824  int ipLo = 0;
825  int ipHi = 0;
826  double fupsilon = 0.;
827  double initTemp = 4.0;
828  double finalTemp = 9.1;
829  double stepTemp = 0.2;
830  fprintf( save.ipPnunit[ipPun],"Species\tLo\tHi\tWLAng");
831  for(double logtemp = initTemp;logtemp < finalTemp;logtemp = logtemp + stepTemp )
832  {
833  fprintf( save.ipPnunit[ipPun],"\t%2.1f",logtemp);
834  }
835  fprintf( save.ipPnunit[ipPun],"\n");
836  for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
837  {
838  for( EmissionList::iterator tr=dBaseTrans[ipSpecies].Emis().begin();
839  tr != dBaseTrans[ipSpecies].Emis().end(); ++tr)
840  {
841  if( dBaseTrans[ipSpecies].chLabel() == "Chianti" )
842  {
843  ipLo = tr->Tran().ipLo();
844  ipHi = tr->Tran().ipHi();
845  fprintf( save.ipPnunit[ipPun],"%s\t%i\t%i\t",
846  dBaseSpecies[ipSpecies].chLabel,ipLo+1,ipHi+1);
847  fprintf( save.ipPnunit[ipPun],"%5.3e",tr->Tran().WLAng());
848  for(double logtemp = initTemp;logtemp < finalTemp;logtemp = logtemp + stepTemp )
849  {
850  fupsilon = CHIANTI_Upsilon(ipSpecies, ipELECTRON, ipHi, ipLo, pow(10.,logtemp));
851  fprintf( save.ipPnunit[ipPun],"\t%5.3e",fupsilon);
852  }
853  fprintf( save.ipPnunit[ipPun],"\n");
854  }
855  }
856  }
857  }
858  }
859 
860  else if( strcmp(save.chSave[ipPun],"COOL") == 0 ||
861  strcmp(save.chSave[ipPun],"EACH") == 0 )
862  {
863  /* save cooling, routine in file of same name */
864  if( strcmp(chTime,"LAST") != 0 )
865  CoolSave(save.ipPnunit[ipPun], save.chSave[ipPun]);
866  }
867 
868  else if( strcmp(save.chSave[ipPun],"DOMI") == 0 )
869  {
870  /* save dominant rates */
871  if( strcmp(chTime,"LAST") != 0 )
872  {
873  molecule *debug_species = findspecies( save.chSpeciesDominantRates[ipPun] );
874  if (debug_species == null_mole)
875  {
876  fprintf( ioQQQ,"Error in SAVE DOMINANT RATES, species %s not found\n",
877  save.chSpeciesDominantRates[ipPun]);
878  }
879  else
880  {
881  fprintf( save.ipPnunit[ipPun],
882  "%e\t%e\t", radius.depth_mid_zone, mole.species[ debug_species->index ].column );
883  mole_dominant_rates( debug_species, save.ipPnunit[ipPun] );
884  }
885  }
886  }
887 
888  else if( strncmp(save.chSave[ipPun],"DYN" , 3) == 0 )
889  {
890  /* save dynamics xxx, information about dynamical solutions */
891  if( strcmp(chTime,"LAST") != 0 )
892  DynaSave(save.ipPnunit[ipPun] ,save.chSave[ipPun][3] );
893  }
894 
895  else if( strcmp(save.chSave[ipPun],"ENTH") == 0 )
896  {
897  if( strcmp(chTime,"LAST") != 0 )
898  fprintf( save.ipPnunit[ipPun],
899  "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
900  radius.depth_mid_zone,
901  phycon.EnthalpyDensity,
902  phycon.EnergyExcitation,
903  phycon.EnergyIonization,
904  phycon.EnergyBinding ,
905  0.5*POW2(wind.windv)*dense.xMassDensity , /* KE */
906  5./2.*pressure.PresGasCurr , /* thermal plus PdV work */
907  magnetic.EnthalpyDensity); /* magnetic terms */
908  }
909 
910  else if( strcmp(save.chSave[ipPun],"COLU") == 0 )
911  {
912  /* save column densities */
913  if( strcmp(chTime,"LAST") == 0 )
914  {
915  PrtColumns(save.ipPnunit[ipPun],"TABLE",ipPun);
916  }
917  }
918 
919  else if( strcmp(save.chSave[ipPun],"COLS") == 0 )
920  {
921  /* save some column densities */
922  if( strcmp(chTime,"LAST") == 0 )
923  {
924  /*TODO unify following with PrtColumns above */
925  save_colden( save.ipPnunit[ipPun] );
926  }
927  }
928 
929  else if( strcmp(save.chSave[ipPun],"COMP") == 0 )
930  {
931  /* Compton energy exchange coefficients */
932  if( strcmp(chTime,"LAST") != 0 )
933  {
934  for( jj=0; jj<rfield.nflux; jj = jj + save.ncSaveSkip)
935  {
936  fprintf( save.ipPnunit[ipPun], "%10.2e%10.2e%10.2e\n",
937  rfield.anu[jj], rfield.comup[jj]/rfield.widflx[jj],
938  rfield.comdn[jj]/rfield.widflx[jj] );
939  }
940  }
941  }
942 
943  /* save continuum commands */
944  else if( strcmp(save.chSave[ipPun],"CON ") == 0 )
945  {
946  /* this is the usual "save continuum" case */
947  /* >>chng 06 apr 03, add every option to do every zone */
948  /* if save every is set then nSaveEveryZone must be positive
949  * was init to -1 */
950  bool lgPrintThis =false;
951  if( save.lgSaveEveryZone[ipPun] )
952  {
953  /* this branch, every option is on line so want to print every n zone */
954  if( strcmp(chTime,"LAST") != 0 )
955  {
956  /* not last zone - print first and intermediate cases */
957  if( nzone==1 )
958  {
959  lgPrintThis = true;
960  }
961  else if( nzone%save.nSaveEveryZone[ipPun]==0 )
962  {
963  lgPrintThis = true;
964  }
965  }
966  else
967  {
968  /* this is last zone, print only if did not trip on above */
969  if( nzone!=1 && nzone%save.nSaveEveryZone[ipPun]!=0 )
970  {
971  lgPrintThis = true;
972  }
973  }
974  }
975  else
976  {
977  /* this branch, not "every", so only print the last zone */
978  if( strcmp(chTime,"LAST") == 0 )
979  lgPrintThis = true;
980  }
981  ASSERT( !save.lgSaveEveryZone[ipPun] || save.nSaveEveryZone[ipPun]>0 );
982  if( lgPrintThis )
983  {
984  if( save.lgSaveEveryZone[ipPun] && nzone!=1)
985  fprintf( save.ipPnunit[ipPun], "%s\n",
986  save.chHashString );
987 
988  /* option to also print same arrays but for cumulative arrays */
989  int nEmType = (int)save.punarg[ipPun][0];
990  ASSERT( nEmType==0 || nEmType==1 );
991 
992  const realnum *trans_coef_total=rfield.getCoarseTransCoef();
993  for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
994  {
995  /* four continua predicted here;
996  * incident, attenuated incident, emitted,
997  * then attenuated incident + emitted, last reflected
998  * reflected continuum is stored relative to inner radius
999  * others are stored for this radius */
1000 
1001  /* NB this code also used in save emitted,
1002  * transmitted continuum commands */
1003 
1004  /* the incident continuum */
1005  flxin = ((double)rfield.flux_total_incident[nEmType][j])
1006  *rfield.anu2[j]*
1007  EN1RYD/rfield.widflx[j];
1008 
1009  // a value < 0. indicates that energy should be conserved
1010  realnum resolution = ( save.Resolution < realnum(0.) ) ?
1011  rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1012 
1013  /* the reflected continuum */
1014  flxref = (rfield.anu2[j]*((double)rfield.ConRefIncid[nEmType][j]+rfield.ConEmitReflec[nEmType][j])/rfield.widflx[j] +
1015  rfield.anu[j]*resolution*rfield.reflin[nEmType][j])*EN1RYD;
1016 
1017  /* the attenuated incident continuum */
1018  flxatt = ((double)rfield.flux[nEmType][j])*rfield.anu2[j]*EN1RYD/
1019  rfield.widflx[j]*radius.r1r0sq *
1020  PrettyTransmission( j, trans_coef_total[j] );
1021 
1022  /* the outward emitted continuum */
1023  conem = ((double)rfield.ConEmitOut[nEmType][j]/
1024  rfield.widflx[j]*rfield.anu2[j] + resolution*
1025  rfield.outlin[nEmType][j]*rfield.anu[j])*radius.r1r0sq*
1026  EN1RYD*geometry.covgeo;
1027 
1028  /* sum of emitted and transmitted continua */
1029  flxtrn = conem + flxatt;
1030 
1031  /* photon energy in appropriate energy or wavelength units */
1032  fprintf( save.ipPnunit[ipPun],"%.5e\t", AnuUnit(rfield.AnuOrg[j]) );
1033  /* incident continuum */
1034  fprintf( save.ipPnunit[ipPun],"%.3e\t", flxin );
1035  /* trans cont */
1036  fprintf( save.ipPnunit[ipPun],"%.3e\t", flxatt );
1037  /* DiffOut cont */
1038  fprintf( save.ipPnunit[ipPun],"%.3e\t", conem );
1039  /* net trans cont */
1040  fprintf( save.ipPnunit[ipPun],"%.3e\t", flxtrn );
1041  /* reflected cont */
1042  fprintf( save.ipPnunit[ipPun],"%.3e\t", flxref );
1043  /* total cont */
1044  fprintf( save.ipPnunit[ipPun],"%.3e\t", flxref + flxtrn );
1045 
1046  /* reflected lines */
1047  fprintf( save.ipPnunit[ipPun],"%.3e\t", rfield.anu[j]*resolution*rfield.reflin[nEmType][j]*EN1RYD );
1048 
1049  /* outward lines */
1050  fprintf( save.ipPnunit[ipPun],"%.3e\t", rfield.anu[j]*resolution*rfield.outlin[nEmType][j]*EN1RYD*
1051  radius.r1r0sq*geometry.covgeo );
1052 
1053  fprintf( save.ipPnunit[ipPun], "%4.4s\t%4.4s\t",
1054  /* line label */
1055  rfield.chLineLabel[j] ,
1056  /* cont label*/
1057  rfield.chContLabel[j] );
1058  /* number of lines within that cell over cell width
1059  * save raw continuum has number by itself */
1060  fprintf( save.ipPnunit[ipPun], "%.2f\n", rfield.line_count[j]/rfield.widflx[j]*rfield.anu[j] );
1061  }
1062  }
1063  }
1064 
1065  /* this is the save spectrum command,
1066  * the new continuum command that will replace the previous one */
1067  else if( strcmp(save.chSave[ipPun],"CONN") == 0 )
1068  {
1069  if( strcmp(chTime,"LAST") == 0 )
1070  SaveNewContinuum( save.ipPnunit[ipPun]);
1071  }
1072 
1073  else if( strcmp(save.chSave[ipPun],"CONC") == 0 )
1074  {
1075  /* save incident continuum */
1076  /* set pointer for possible change in units of energy in continuum
1077  * AnuUnit will give anu in whatever units were set with save units */
1078  if( strcmp(chTime,"LAST") == 0 )
1079  {
1080  /* incident continuum */
1081  for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
1082  {
1083  flxin = rfield.flux_total_incident[0][j]*rfield.anu2[j]*
1084  EN1RYD/rfield.widflx[j];
1085  /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */
1086  fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\n",
1087  AnuUnit(rfield.AnuOrg[j]), flxin );
1088  }
1089  }
1090  }
1091 
1092  else if( strcmp(save.chSave[ipPun],"CONG") == 0 )
1093  {
1094  /* save emitted grain continuum in optically thin limit */
1095  if( strcmp(chTime,"LAST") == 0 )
1096  {
1097  /* GrainMakeDiffuse broke out emission into types
1098  * according to matType */
1099  for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
1100  {
1101  fsum = gv.GraphiteEmission[j]*rfield.anu2[j]*
1102  EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo;
1103 
1104  fout = gv.SilicateEmission[j]*rfield.anu2[j]*
1105  EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo;
1106 
1107  /* anu is .5e format to resolve energy mesh near 1 Ryd
1108  * AnuUnit gives anu in whatever units were set with units option */
1109  fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\n",
1110  AnuUnit(rfield.AnuOrg[j]) , fsum , fout ,
1111  /* total emission per unit volume defined in GrainMakeDiffuse
1112  * used in RT_diffuse to form total grain emission */
1113  gv.GrainEmission[j]*rfield.anu2[j]*
1114  EN1RYD/rfield.widflx[j] *radius.r1r0sq*geometry.covgeo );
1115  }
1116  }
1117  }
1118 
1119  else if( strcmp(save.chSave[ipPun],"CONR") == 0 )
1120  {
1121  /* save reflected continuum */
1122  /* set pointer for possible change in units of energy in continuum
1123  * AnuUnit will give anu in whatever units were set with save units */
1124  if( strcmp(chTime,"LAST") == 0 )
1125  {
1126  if( geometry.lgSphere )
1127  {
1128  fprintf( save.ipPnunit[ipPun], " Reflected continuum not predicted when SPHERE is set.\n" );
1129  fprintf( ioQQQ ,
1130  "\n\n>>>>>>>>>>>>>\n Reflected continuum not predicted when SPHERE is set.\n" );
1131  cdEXIT(EXIT_FAILURE);
1132  }
1133 
1134  for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
1135  {
1136  // a value < 0. indicates that energy should be conserved
1137  realnum resolution = ( save.Resolution < realnum(0.) ) ?
1138  rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1139 
1140  /* reflected continuum */
1141  flxref = rfield.anu2[j]*((double)rfield.ConRefIncid[0][j]+rfield.ConEmitReflec[0][j])/
1142  rfield.widflx[j]*EN1RYD;
1143  /* reflected lines */
1144  fref = rfield.anu[j]*resolution*rfield.reflin[0][j]*EN1RYD;
1145  /* ratio of reflected to incident continuum, the albedo */
1146  if( rfield.flux_total_incident[0][j] > 1e-25 )
1147  {
1148  av = rfield.ConRefIncid[0][j]/rfield.flux_total_incident[0][j];
1149  }
1150  else
1151  {
1152  av = 0.;
1153  }
1154  /* >>chng 96 oct 22, format of anu to .5 to resolve energy mesh near 1 Ryd */
1155  fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4s\n",
1156  AnuUnit(rfield.AnuOrg[j]), flxref, fref, flxref + fref,
1157  av, rfield.chContLabel[j] );
1158  }
1159  }
1160  }
1161 
1162  else if( strcmp(save.chSave[ipPun],"CNVE") == 0 )
1163  {
1164  /* the save convergence error command */
1165  if( strcmp(chTime,"LAST") != 0 )
1166  {
1167  fprintf( save.ipPnunit[ipPun],
1168  "%.5e\t%li\t%.4e\t%.4f\t%.4e\t%.4e\t%.3f\t%.4e\t%.4e\t%.4f\n",
1169  radius.depth_mid_zone,
1170  conv.nPres2Ioniz,
1171  pressure.PresTotlCurr,
1172  pressure.PresTotlError*100.,
1173  dense.EdenTrue,
1174  dense.eden,
1175  (dense.EdenTrue - dense.eden)*100./dense.EdenTrue,
1176  thermal.htot,
1177  thermal.ctot,
1178  (thermal.htot - thermal.ctot)*100./thermal.htot );
1179  }
1180  }
1181 
1182  else if( strcmp(save.chSave[ipPun],"CONB") == 0 )
1183  {
1184  /* save continuum bins binning */
1185  /* set pointer for possible change in units of energy in continuum
1186  * AnuUnit will give anu in whatever units were set with save units */
1187  if( strcmp(chTime,"LAST") != 0 )
1188  {
1189  for( j=0; j<rfield.nupper; j = j + save.ncSaveSkip)
1190  {
1191  fprintf( save.ipPnunit[ipPun], "%14.5e %14.5e %14.5e\n",
1192  AnuUnit(rfield.AnuOrg[j]), rfield.anu[j], rfield.widflx[j] );
1193  }
1194  }
1195  }
1196 
1197  else if( strcmp(save.chSave[ipPun],"COND") == 0 )
1198  {
1199  ASSERT( fp_equal( save.punarg[ipPun][0] , (realnum)2. ) ||
1200  fp_equal( save.punarg[ipPun][0] , (realnum)1.) );
1201 
1202  /* save diffuse continuum the local line and continuous emission */
1203  if( strcmp(chTime,"LAST") == 0 &&
1204  fp_equal(save.punarg[ipPun][0] , (realnum)1.) )
1205  {
1206  /* this option to save diffuse emission for last zone
1207  * as series of columns */
1208  for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
1209  {
1210  // a value < 0. indicates that energy should be conserved
1211  realnum resolution = ( save.Resolution < realnum(0.) ) ?
1212  rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1213 
1214  double EmisLin = resolution*EN1RYD*
1215  rfield.DiffuseLineEmission[j]*rfield.anu[j];
1216  double EmisCon = rfield.ConEmitLocal[nzone][j]*
1217  rfield.anu2[j]*EN1RYD/rfield.widflx[j];
1218  fprintf( save.ipPnunit[ipPun], "%.5e\t%.5e\t%.5e\t%.5e\n",
1219  AnuUnit(rfield.AnuOrg[j]),
1220  EmisCon ,
1221  EmisLin ,
1222  EmisCon+EmisLin);
1223  }
1224  }
1225  else if( strcmp(chTime,"LAST") != 0 &&
1226  fp_equal(save.punarg[ipPun][0] , (realnum)2.) )
1227  {
1228  /* diffuse emission per zone, with each a long row */
1229  static bool lgMustPrintHeader=true;
1230  if( lgMustPrintHeader )
1231  {
1232  lgMustPrintHeader = false;
1233  fprintf( save.ipPnunit[ipPun], "%.5e",
1234  AnuUnit(rfield.AnuOrg[0]) );
1235  for( j=1; j<rfield.nflux; j = j+save.ncSaveSkip)
1236  {
1237  fprintf( save.ipPnunit[ipPun], "\t%.5e",
1238  AnuUnit(rfield.AnuOrg[j]) );
1239  }
1240  fprintf( save.ipPnunit[ipPun], "\n" );
1241  }
1242  // a value < 0. indicates that energy should be conserved
1243  realnum resolution = ( save.Resolution < realnum(0.) ) ?
1244  rfield.anu[0]/rfield.widflx[0] : save.Resolution;
1245  double EmisLin = resolution*EN1RYD*
1246  rfield.DiffuseLineEmission[0]*rfield.anu[0];
1247  double EmisCon = rfield.ConEmitLocal[nzone][0]*
1248  rfield.anu2[0]*EN1RYD/rfield.widflx[0];
1249  fprintf( save.ipPnunit[ipPun], "%.5e",
1250  EmisCon+EmisLin);
1251  for( j=1; j<rfield.nflux; j = j+save.ncSaveSkip)
1252  {
1253  // a value < 0. indicates that energy should be conserved
1254  resolution = ( save.Resolution < realnum(0.) ) ?
1255  rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1256  double EmisLin = resolution*EN1RYD*
1257  rfield.DiffuseLineEmission[j]*rfield.anu[j];
1258  double EmisCon = rfield.ConEmitLocal[nzone][j]*
1259  rfield.anu2[j]*EN1RYD/rfield.widflx[j];
1260  fprintf( save.ipPnunit[ipPun], "\t%.5e",
1261  EmisCon+EmisLin);
1262  }
1263  fprintf( save.ipPnunit[ipPun], "\n" );
1264  }
1265  }
1266 
1267  else if( strcmp(save.chSave[ipPun],"CONE") == 0 )
1268  {
1269  /* save emitted continuum */
1270  /* set pointer for possible change in units of energy in continuum
1271  * AnuUnit will give anu in whatever units were set with save units */
1272  if( strcmp(chTime,"LAST") == 0 )
1273  {
1274  /* save emitted continuum */
1275  for( j=0; j<rfield.nflux;j = j +save.ncSaveSkip)
1276  {
1277  // a value < 0. indicates that energy should be conserved
1278  realnum resolution = ( save.Resolution < realnum(0.) ) ?
1279  rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1280 
1281  /* this is the reflected component */
1282  flxref = (rfield.anu2[j]*(rfield.ConRefIncid[0][j]+rfield.ConEmitReflec[0][j])/
1283  rfield.widflx[j] + rfield.anu[j]*resolution*
1284  rfield.reflin[0][j])*EN1RYD;
1285 
1286  /* this is the total emission in the outward direction */
1287  conem = (rfield.ConEmitOut[0][j])/rfield.widflx[j]*rfield.anu2[j] +
1288  resolution*rfield.outlin[0][j]*rfield.anu[j];
1289 
1290  conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;
1291 
1292  /* output: photon energy, reflected, outward, total emission
1293  * >>chng 96 oct 22, format of anu to .5e to resolve energy mesh near 1 Ryd */
1294  fprintf( save.ipPnunit[ipPun], "%.5e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\n",
1295  AnuUnit(rfield.AnuOrg[j]),
1296  flxref,
1297  conem,
1298  flxref + conem,
1299  rfield.chLineLabel[j],
1300  rfield.chContLabel[j]
1301  );
1302  }
1303  }
1304  }
1305 
1306  /* save fine continuum command */
1307  else if( strcmp(save.chSave[ipPun],"CONf") == 0 )
1308  {
1309  if( strcmp(chTime,"LAST") == 0 )
1310  {
1311  long nu_hi , nskip;
1312  if( save.punarg[ipPun][0] > 0. )
1313  /* possible lower bounds to energy range -
1314  * 0 if not set with range option*/
1315  j = ipFineCont( save.punarg[ipPun][0] );
1316  else
1317  j = 0;
1318 
1319  /* upper limit set with range option */
1320  if( save.punarg[ipPun][1]> 0. )
1321  nu_hi = ipFineCont( save.punarg[ipPun][1]);
1322  else
1323  nu_hi = rfield.nfine;
1324 
1325  /* number of cells to bring together, default is 10 */
1326  nskip = (long)save.punarg[ipPun][2];
1327 
1328  do
1329  {
1330  realnum sum1 = rfield.fine_opt_depth[j];
1331  realnum xnu = rfield.fine_anu[j];
1332  for( jj=1; jj<nskip; ++jj )
1333  {
1334  xnu += rfield.fine_anu[j+jj];
1335  sum1 += rfield.fine_opt_depth[j+jj];
1336  }
1337  fprintf( save.ipPnunit[ipPun],
1338  "%.6e\t%.3e\n",
1339  AnuUnit(xnu/nskip),
1340  sexp(sum1/nskip) );
1341  j += nskip;
1342  } while( j < nu_hi );
1343  }
1344  }
1345 
1346  else if( strcmp(save.chSave[ipPun],"CONi") == 0 )
1347  {
1348  /* save continuum interactions */
1349  /* set pointer for possible change in units of energy in continuum
1350  * AnuUnit will give anu in whatever units were set with save units */
1351 
1352  /* continuum interactions */
1353  if( strcmp(chTime,"LAST") != 0 )
1354  {
1355  /* this is option to set lowest energy */
1356  if( save.punarg[ipPun][0] <= 0. )
1357  {
1358  i1 = 1;
1359  }
1360  else if( save.punarg[ipPun][0] < 100. )
1361  {
1362  i1 = ipoint(save.punarg[ipPun][0]);
1363  }
1364  else
1365  {
1366  i1 = (long int)save.punarg[ipPun][0];
1367  }
1368 
1369  fref = 0.;
1370  fout = 0.;
1371  fsum = 0.;
1372  sum = 0.;
1373  flxin = 0.;
1374 
1375  for( j=i1-1; j < rfield.nflux; j++ )
1376  {
1377  fref += rfield.flux[0][j]*opac.opacity_abs[j];
1378  fout += rfield.otslin[j]*opac.opacity_abs[j];
1379  fsum += rfield.otscon[j]*opac.opacity_abs[j];
1380  sum += rfield.ConInterOut[j]*opac.opacity_abs[j];
1381  flxin += (rfield.outlin[0][j] + rfield.outlin_noplot[j])*opac.opacity_abs[j];
1382  }
1383  fprintf( save.ipPnunit[ipPun], "%10.2e%10.2e%10.2e%10.2e%10.2e\n",
1384  fref, fout, fsum, sum, flxin );
1385  }
1386  }
1387 
1388  else if( strcmp(save.chSave[ipPun],"CONI") == 0 )
1389  {
1390  /* save ionizing continuum */
1391  /* set pointer for possible change in units of energy in continuum
1392  * AnuUnit will give anu in whatever units were set with save units */
1393 
1394  if( save.lgSaveEveryZone[ipPun] || (strcmp(chTime,"LAST") == 0) )
1395  {
1396  /* this flag will remember whether we have ever printed anything */
1397  bool lgPrt=false;
1398  if( save.lgSaveEveryZone[ipPun] )
1399  fprintf(save.ipPnunit[ipPun],"#save every zone %li\n", nzone);
1400 
1401  /* save ionizing continuum command
1402  * this is option to set lowest energy,
1403  * if no number was entered then this was zero */
1404  if( save.punarg[ipPun][0] <= 0. )
1405  i1 = 1;
1406  else if( save.punarg[ipPun][0] < 100. )
1407  i1 = ipoint(save.punarg[ipPun][0]);
1408  else
1409  i1 = (long int)save.punarg[ipPun][0];
1410 
1411  sum = 0.;
1412  for( j=i1-1; j < rfield.nflux; j++ )
1413  {
1414  flxcor = rfield.flux[0][j] +
1415  rfield.otslin[j] +
1416  rfield.otscon[j] +
1417  rfield.ConInterOut[j] +
1418  rfield.outlin[0][j] + rfield.outlin_noplot[j];
1419 
1420  sum += flxcor*opac.opacity_abs[j];
1421  }
1422 
1423  if( sum > 0. )
1424  sum = 1./sum;
1425  else
1426  sum = 1.;
1427 
1428  fsum = 0.;
1429 
1430  for( j=i1-1; j<rfield.nflux; ++j)
1431  {
1432  flxcor = rfield.flux[0][j] +
1433  rfield.otslin[j] +
1434  rfield.otscon[j] +
1435  rfield.ConInterOut[j]+
1436  rfield.outlin[0][j] + rfield.outlin_noplot[j];
1437 
1438  fsum += flxcor*opac.opacity_abs[j];
1439 
1440  /* punched quantities are freq, flux, flux*cross sec,
1441  * fraction of total, integral fraction of total */
1442  RateInter = flxcor*opac.opacity_abs[j]*sum;
1443 
1444  /* punage(ipPun,2) is lowest interaction rate to consider, def=0.01 (1 percent) */
1445  /* >>chng 01 nov 22, format to c-friendly */
1446  if( (RateInter >= save.punarg[ipPun][1]) && (flxcor > SMALLFLOAT) )
1447  {
1448  lgPrt = true;
1449  /* >>chng 96 oct 22, format of anu to 11.5 to resolve energy mesh near 1 Ryd */
1450  fprintf( save.ipPnunit[ipPun],
1451  "%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",
1452  j,
1453  AnuUnit(rfield.AnuOrg[j]),
1454  flxcor,
1455  flxcor*opac.opacity_abs[j],
1456  rfield.flux[0][j]/flxcor,
1457  rfield.otslin[j]/flxcor,
1458  rfield.otscon[j]/flxcor,
1459  (rfield.outlin[0][j] + rfield.outlin_noplot[j])/flxcor,
1460  rfield.ConInterOut[j]/flxcor,
1461  RateInter,
1462  fsum*sum,
1463  rfield.chLineLabel[j],
1464  rfield.chContLabel[j] );
1465  }
1466  }
1467  if( !lgPrt )
1468  {
1469  /* entered logical block but did not print anything */
1470  fprintf(save.ipPnunit[ipPun],
1471  " punchdo, the PUNCH IONIZING CONTINUUM command "
1472  "did not find a strong point, sum and fsum were %.2e %.2e\n",
1473  sum,fsum);
1474  fprintf(save.ipPnunit[ipPun],
1475  " punchdo, the low-frequency energy was %.5e Ryd\n",
1476  rfield.anu[i1-1]);
1477  fprintf(save.ipPnunit[ipPun],
1478  " You can reset the threshold for the lowest fractional "
1479  "interaction to print with the second number of the save command\n"
1480  " The fraction was %.3f and this was too large.\n",
1481  save.punarg[ipPun][1]);
1482  }
1483  }
1484  }
1485 #ifdef USE_NLTE7
1486  else if( strcmp(save.chSave[ipPun],"CONl") == 0 )
1487  {
1488  if( strcmp(chTime,"LAST") != 0 )
1489  {
1490  int nEmType = (int)save.punarg[ipPun][0];
1491  ASSERT( nEmType==0 || nEmType==1 );
1492 
1493  for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
1494  {
1495 
1496  // a value < 0. indicates that energy should be conserved
1497  realnum resolution = ( save.Resolution < realnum(0.) ) ?
1498  rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1499 
1500  /* the reflected continuum */
1501  flxref = (rfield.anu2[j]*((double)rfield.ConRefIncid[nEmType][j]+rfield.ConEmitReflec[nEmType][j])/rfield.widflx[j] +
1502  rfield.anu[j]*resolution*rfield.reflin[nEmType][j])*EN1RYD;
1503 
1504  /* the attenuated incident continuum */
1505  flxatt = rfield.flux[nEmType][j]*rfield.anu2[j]*EN1RYD/
1506  rfield.widflx[j]*radius.r1r0sq *
1507  PrettyTransmission( j, rfield.trans_coef_total[j] );
1508 
1509  /* the outward emitted continuum */
1510  conem = (((double)rfield.ConEmitOut[nEmType][j])/
1511  rfield.widflx[j]*rfield.anu2[j] + resolution*
1512  rfield.outlin[nEmType][j]*rfield.anu[j])*radius.r1r0sq*
1513  EN1RYD*geometry.covgeo;
1514 
1515  /* sum of emitted and transmitted continua */
1516  flxtrn = conem + flxatt;
1517 
1518  //Set upper and lower limits on which wavelength/energy values are printed.
1519  double lowlim, highlim, NRGeV;
1520  lowlim = 1.0;
1521  highlim = 250.0;
1522  NRGeV = AnuUnit(rfield.AnuOrg[j]);
1523 
1524  if( NRGeV >= lowlim && NRGeV <= highlim )
1525  {
1526  /* photon energy in appropriate energy or wavelength units */
1527  fprintf( save.ipPnunit[ipPun],"%14.8e ", NRGeV );
1528  /* print zeroes for BB BF and FF */
1529  fprintf( save.ipPnunit[ipPun],"%14.8e %14.8e %14.8e ",0.,0.,0.);
1530  /* total cont */
1531  fprintf( save.ipPnunit[ipPun],"%14.8e\n", (flxref + flxtrn)/NRGeV );
1532  }
1533  }
1534  }
1535 
1536  }
1537 #endif
1538 
1539 
1540  else if( strcmp(save.chSave[ipPun],"CONS") == 0 )
1541  {
1542  if( strcmp(chTime,"LAST") != 0 )
1543  {
1544  // continuum volume emissivity and opacity as a function of radius
1545  // command was "save continuum emissivity" */
1546  if( save.ipEmisFreq[ipPun] < 0 )
1547  save.ipEmisFreq[ipPun] = ipoint(save.emisfreq[ipPun].Ryd());
1548  j = save.ipEmisFreq[ipPun]-1;
1549 
1550  fprintf( save.ipPnunit[ipPun],
1551  "%.14e\t%.14e\t%.5e\t%.5e\t%.5e\n",
1552  radius.Radius_mid_zone,
1553  radius.depth_mid_zone,
1554  rfield.anu2[j]*rfield.ConEmitLocal[nzone][j]/rfield.widflx[j]*EN1RYD,
1555  opac.opacity_abs[j],
1556  opac.opacity_sct[j] );
1557  }
1558  }
1559 
1560  else if( strcmp(save.chSave[ipPun],"CORA") == 0 )
1561  {
1562  /* save raw continuum */
1563  /* set pointer for possible change in units of energy in continuum
1564  * AnuUnit will give anu in whatever units were set with save units */
1565 
1566  if( strcmp(chTime,"LAST") == 0 )
1567  {
1568  /* this option to save all raw ionizing continuum */
1569  for( j=0;j<rfield.nflux;j = j + save.ncSaveSkip)
1570  {
1571  fprintf( save.ipPnunit[ipPun],
1572  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%4.4s\t%4.4s\t",
1573  AnuUnit(rfield.AnuOrg[j]),
1574  rfield.flux[0][j],
1575  rfield.otslin[j],
1576  rfield.otscon[j],
1577  rfield.ConRefIncid[0][j],
1578  rfield.ConEmitReflec[0][j],
1579  rfield.ConInterOut[j],
1580  rfield.outlin[0][j]+rfield.outlin_noplot[j],
1581  rfield.ConEmitOut[0][j],
1582  rfield.chLineLabel[j],
1583  rfield.chContLabel[j]
1584  );
1585  /* number of lines within that cell */
1586  fprintf( save.ipPnunit[ipPun], "%li\n", rfield.line_count[j] );
1587  }
1588  }
1589  }
1590 
1591  else if( strcmp(save.chSave[ipPun],"CONo") == 0 )
1592  {
1593  /* save outward local continuum */
1594  /* set pointer for possible change in units of energy in continuum
1595  * AnuUnit will give anu in whatever units were set with save units */
1596 
1597  if( strcmp(chTime,"LAST") == 0 )
1598  {
1599  /* option to save out outward continuum here */
1600  for( j=0;j<rfield.nflux; j = j + save.ncSaveSkip)
1601  {
1602  fprintf( save.ipPnunit[ipPun], "%11.5e%10.2e%10.2e\n",
1603  AnuUnit(rfield.AnuOrg[j]),
1604  rfield.ConEmitOut[0][j]+ rfield.outlin[0][j] + rfield.outlin_noplot[j],
1605  (rfield.flux[0][j] + rfield.otscon[j] + rfield.otslin[j] +
1606  rfield.ConInterOut[j])*opac.opacity_abs[j]*
1607  rfield.anu[j] );
1608  }
1609  }
1610  }
1611 
1612  else if( strcmp(save.chSave[ipPun],"CONO") == 0 )
1613  {
1614  /* save outward continuum */
1615  /* set pointer for possible change in units of energy in continuum
1616  * AnuUnit will give anu in whatever units were set with save units */
1617 
1618  if( strcmp(chTime,"LAST") == 0 )
1619  {
1620  /* option to save out outward continuum */
1621  for( j=0; j<rfield.nflux; j = j + save.ncSaveSkip)
1622  {
1623  // a value < 0. indicates that energy should be conserved
1624  realnum resolution = ( save.Resolution < realnum(0.) ) ?
1625  rfield.anu[j]/rfield.widflx[j] : save.Resolution;
1626 
1627  fprintf( save.ipPnunit[ipPun], "%11.5e%10.2e%10.2e%10.2e%10.2e\n",
1628  AnuUnit(rfield.AnuOrg[j]),
1629  rfield.flux[0][j]*rfield.anu2[j]* EN1RYD/rfield.widflx[j]*radius.r1r0sq,
1630  rfield.ConInterOut[j]/rfield.widflx[j]*rfield.anu2[j]*radius.r1r0sq*EN1RYD,
1631  resolution*(rfield.outlin[0][j]+rfield.outlin_noplot[j])*rfield.anu[j]*radius.r1r0sq*EN1RYD,
1632  (rfield.ConInterOut[j]/rfield.widflx[j]*
1633  rfield.anu2[j] + resolution*(rfield.outlin[0][j]+rfield.outlin_noplot[j])*
1634  rfield.anu[j])*radius.r1r0sq*EN1RYD );
1635  }
1636  }
1637  }
1638 
1639  else if( strcmp(save.chSave[ipPun],"CONT") == 0 )
1640  {
1641  /* save transmitted continuum - this is not the main "save continuum"
1642  * command - search on "CON " above
1643  * set pointer for possible change in units of energy in continuum
1644  * AnuUnit will give anu in whatever units were set with save units */
1645 
1646  if( strcmp(chTime,"LAST") == 0 )
1647  {
1648  fprintf( save.ipPnunit[ipPun], "#\n" );
1649  fprintf( save.ipPnunit[ipPun], "%32ld # file format version number\n",
1650  VERSION_TRNCON );
1651  fprintf( save.ipPnunit[ipPun], "%s # check 1\n",
1652  continuum.mesh_md5sum.c_str() );
1653  union {
1654  double x;
1655  uint32 i[2];
1656  } u;
1657  u.x = double(rfield.emm);
1658  if( cpu.i().big_endian() )
1659  fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 2\n",
1660  u.i[0], u.i[1] );
1661  else
1662  fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 2\n",
1663  u.i[1], u.i[0] );
1664  u.x = double(rfield.egamry);
1665  if( cpu.i().big_endian() )
1666  fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 3\n",
1667  u.i[0], u.i[1] );
1668  else
1669  fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 3\n",
1670  u.i[1], u.i[0] );
1671  u.x = continuum.ResolutionScaleFactor;
1672  if( cpu.i().big_endian() )
1673  fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 4\n",
1674  u.i[0], u.i[1] );
1675  else
1676  fprintf( save.ipPnunit[ipPun], "%23.8x %8.8x # check 4\n",
1677  u.i[1], u.i[0] );
1678  fprintf( save.ipPnunit[ipPun], "%32ld # nflux\n",
1679  (rfield.nflux+save.ncSaveSkip-1)/save.ncSaveSkip );
1680  fprintf( save.ipPnunit[ipPun], "#\n" );
1681 
1682  const realnum *trans_coef_total=rfield.getCoarseTransCoef();
1683 
1684  /* this option to save transmitted continuum */
1685  for( j=0; j < rfield.nflux; j += save.ncSaveSkip )
1686  {
1687  /* attenuated incident continuum
1688  * >>chng 97 jul 10, remove SaveLWidth from this one only since
1689  * we must conserve energy even in lines
1690  * >>chng 07 apr 26 include transmission coefficient */
1691  flxatt = rfield.flux[0][j]*rfield.anu2[j]*EN1RYD/
1692  rfield.widflx[j]*radius.r1r0sq*trans_coef_total[j];
1693 
1694  /*conem = (rfield.ConOutNoInter[j] + rfield.ConInterOut[j]+rfield.outlin[0][j])*
1695  rfield.anu2[j];
1696  conem *= radius.r1r0sq*EN1RYD*geometry.covgeo;*/
1697  /* >>chng 00 jan 03, above did not include all contributors.
1698  * Pasted in below from usual
1699  * save continuum command */
1700  /* >>chng 04 jul 15, removed factor of save.SaveLWidth -
1701  * this should not be there to conserve energy, as explained in hazy
1702  * where command was documented, and in comment above. caught by PvH */
1703  /* >>chng 04 jul 23, incorrect use of outlin - before multiplied by an2,
1704  * quantity should be photons per Ryd, since init quantity is
1705  * photons per cell. Must div by widflx. caught by PvH */
1706  /*conem = (rfield.ConEmitOut[0][j]/rfield.widflx[j]*rfield.anu2[j] +
1707  rfield.outlin[0][j]*rfield.anu[j])*radius.r1r0sq*
1708  EN1RYD*geometry.covgeo;*/
1709  conem = (rfield.ConEmitOut[0][j] + rfield.outlin[0][j]) / rfield.widflx[j]*
1710  rfield.anu2[j]*radius.r1r0sq*EN1RYD*geometry.covgeo;
1711 
1712  flxtrn = conem + flxatt;
1713 
1714  /* use AnuOrg here instead of anu since probably
1715  * going to be used by table read
1716  * and we want good anu array for sanity check*/
1717  fprintf( save.ipPnunit[ipPun],"%.5e\t%.3e\t%.3e\n",
1718  AnuUnit(rfield.AnuOrg[j]), flxtrn,
1719  trans_coef_total[j] );
1720  }
1721  }
1722  }
1723 
1724  else if( strcmp(save.chSave[ipPun],"CON2") == 0 )
1725  {
1726  /* save total two-pohoton continuum */
1727  if( strcmp(chTime,"LAST") == 0 )
1728  {
1729  /* this option to save diffuse continuum */
1730  for( j=0; j<rfield.nflux; j = j+save.ncSaveSkip)
1731  {
1732  fprintf( save.ipPnunit[ipPun], "%.5e\t%.5e\t%.5e\n",
1733  AnuUnit(rfield.AnuOrg[j]),
1734  rfield.TotDiff2Pht[j]/rfield.widflx[j] ,
1735  rfield.TotDiff2Pht[j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
1736  }
1737  }
1738  }
1739 
1740  else if( strcmp(save.chSave[ipPun],"DUSE") == 0 )
1741  {
1742  /* save grain extinction - includes only grain opacity, not total */
1743  if( strcmp(chTime,"LAST") != 0 )
1744  {
1745  fprintf( save.ipPnunit[ipPun], " %.5e\t",
1746  radius.depth_mid_zone );
1747 
1748  /* visual extinction of an extended source (like a PDR)*/
1749  fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
1750 
1751  /* visual extinction of point source (star)*/
1752  fprintf( save.ipPnunit[ipPun], "%.2e\n" , rfield.extin_mag_V_point);
1753  }
1754  }
1755 
1756  else if( strcmp(save.chSave[ipPun],"DUSO") == 0 )
1757  {
1758  /* save grain opacity */
1759  if( strcmp(chTime,"LAST") == 0 )
1760  {
1761  for( j=0; j < rfield.nflux; j++ )
1762  {
1763  double scat;
1764  fprintf( save.ipPnunit[ipPun],
1765  "%.5e\t%.2e\t%.2e\t%.2e\t",
1766  /* photon energy or wavelength */
1767  AnuUnit(rfield.AnuOrg[j]),
1768  /* total opacity, discount forward scattering */
1769  gv.dstab[j] + gv.dstsc[j],
1770  /* absorption opacity */
1771  gv.dstab[j],
1772  /* scatter, with forward discounted */
1773  gv.dstsc[j] );
1774  /* add together total scattering, discounting 1-g */
1775  scat = 0.;
1776  /* sum over all grain species */
1777  for( size_t nd=0; nd < gv.bin.size(); nd++ )
1778  {
1779  scat += gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->dstAbund;
1780  }
1781  /* finally, scattering including effects of forward scattering */
1782  fprintf( save.ipPnunit[ipPun],
1783  "%.2e\t", scat );
1784  fprintf( save.ipPnunit[ipPun],
1785  "%.2e\n", gv.dstsc[j]/SDIV(gv.dstab[j] + gv.dstsc[j]) );
1786  }
1787  }
1788  }
1789 
1790  /* save grain abundance and save grain D/G ratio commands */
1791  else if( strcmp(save.chSave[ipPun],"DUSA") == 0 ||
1792  strcmp(save.chSave[ipPun],"DUSD") == 0 )
1793  {
1794  bool lgDGRatio = ( strcmp(save.chSave[ipPun],"DUSD") == 0 );
1795 
1796  /* grain abundance */
1797  if( strcmp(chTime,"LAST") != 0 )
1798  {
1799  /* used to print header exactly one time */
1800  static bool lgMustPrtHeaderDRRatio = true,
1801  lgMustPrtHeaderAbundance=true;
1802  /* print grain header first if this has not yet been done */
1803  if( ( lgMustPrtHeaderDRRatio && lgDGRatio ) ||
1804  ( lgMustPrtHeaderAbundance && !lgDGRatio ) )
1805  {
1806  /* only print one header for each case, but must
1807  * track separately if both used in same sim */
1808  if( lgMustPrtHeaderDRRatio && lgDGRatio )
1809  lgMustPrtHeaderDRRatio = false;
1810  else if( lgMustPrtHeaderAbundance &&!lgDGRatio )
1811  lgMustPrtHeaderAbundance = false;
1812 
1813  fprintf( save.ipPnunit[ipPun], "#Depth" );
1814  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1815  fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1816  fprintf( save.ipPnunit[ipPun], "\ttotal\n" );
1817 
1818  /* now print grain radius converting from cm to microns */
1819  fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
1820  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1821  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1822  fprintf( save.ipPnunit[ipPun], "\txxxx\n" );
1823  }
1824  fprintf( save.ipPnunit[ipPun], " %.5e",
1825  radius.depth_mid_zone );
1826  /* grain abundance per bin in g/cm^3 */
1827  double total = 0.;
1828  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1829  {
1830  double abund = gv.bin[nd]->IntVol*gv.bin[nd]->dustp[0]*
1831  gv.bin[nd]->cnv_H_pCM3;
1832  if( lgDGRatio )
1833  abund /= dense.xMassDensity;
1834  fprintf( save.ipPnunit[ipPun], "\t%.3e", abund );
1835  total += abund;
1836  }
1837  fprintf( save.ipPnunit[ipPun], "\t%.3e\n", total );
1838  }
1839  }
1840 
1841  else if( strcmp(save.chSave[ipPun],"DUSP") == 0 )
1842  {
1843  /* grain potential */
1844  if( strcmp(chTime,"LAST") != 0 )
1845  {
1846  /* used to print header exactly one time */
1847  static bool lgMustPrtHeader = true;
1848  /* do labels first if this is first zone */
1849  if( lgMustPrtHeader )
1850  {
1851  /* first print string giving grain id */
1852  fprintf( save.ipPnunit[ipPun], "#Depth" );
1853  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1854  fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1855  fprintf( save.ipPnunit[ipPun], "\n" );
1856 
1857  /* now print grain radius converting from cm to microns */
1858  fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
1859  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1860  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1861  fprintf( save.ipPnunit[ipPun], "\n" );
1862 
1863  /* don't need to do this, ever again */
1864  lgMustPrtHeader = false;
1865  }
1866  fprintf( save.ipPnunit[ipPun], " %.5e",
1867  radius.depth_mid_zone );
1868  /* grain potential in eV */
1869  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1870  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->dstpot*EVRYD );
1871  fprintf( save.ipPnunit[ipPun], "\n" );
1872  }
1873  }
1874 
1875  else if( strcmp(save.chSave[ipPun],"DUSR") == 0 )
1876  {
1877  /* grain H2 formation rates */
1878  if( strcmp(chTime,"LAST") != 0 )
1879  {
1880  /* used to print header exactly one time */
1881  static bool lgMustPrtHeader = true;
1882  /* do labels first if this is first zone */
1883  if( lgMustPrtHeader )
1884  {
1885  /* first print string giving grain id */
1886  fprintf( save.ipPnunit[ipPun], "#Depth" );
1887  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1888  fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1889  fprintf( save.ipPnunit[ipPun], "\n" );
1890 
1891  /* now print grain radius converting from cm to microns */
1892  fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
1893  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1894  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1895  fprintf( save.ipPnunit[ipPun], "\n" );
1896 
1897  /* don't need to do this, ever again */
1898  lgMustPrtHeader = false;
1899  }
1900  fprintf( save.ipPnunit[ipPun], " %.5e",
1901  radius.depth_mid_zone );
1902  /* grain formation rate for H2 */
1903  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1904  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->rate_h2_form_grains_used );
1905  fprintf( save.ipPnunit[ipPun], "\n" );
1906  }
1907  }
1908 
1909  else if( strcmp(save.chSave[ipPun],"DUST") == 0 )
1910  {
1911  /* grain temperatures - K*/
1912  if( strcmp(chTime,"LAST") != 0 )
1913  {
1914  /* used to print header exactly one time */
1915  static bool lgMustPrtHeader = true;
1916  /* do labels first if this is first zone */
1917  if( lgMustPrtHeader )
1918  {
1919  /* first print string giving grain id */
1920  fprintf( save.ipPnunit[ipPun], "#Depth" );
1921  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1922  fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1923  fprintf( save.ipPnunit[ipPun], "\n" );
1924 
1925  /* now print grain radius converting from cm to microns */
1926  fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
1927  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1928  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1929  fprintf( save.ipPnunit[ipPun], "\n" );
1930 
1931  /* don't need to do this, ever again */
1932  lgMustPrtHeader = false;
1933  }
1934  fprintf( save.ipPnunit[ipPun], " %.5e",
1935  radius.depth_mid_zone );
1936  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1937  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->tedust );
1938  fprintf( save.ipPnunit[ipPun], "\n" );
1939  }
1940  }
1941 
1942  else if( strcmp(save.chSave[ipPun],"DUSC") == 0 )
1943  {
1944  /* save grain charge - eden from grains and
1945  * charge per grain in electrons / grain */
1946  if( strcmp(chTime,"LAST") != 0 )
1947  {
1948  /* used to print header exactly one time */
1949  static bool lgMustPrtHeader = true;
1950  /* do labels first if this is first zone */
1951  if( lgMustPrtHeader )
1952  {
1953  /* first print string giving grain id */
1954  fprintf( save.ipPnunit[ipPun], "#Depth\tne(grn)" );
1955  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1956  fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1957  fprintf( save.ipPnunit[ipPun], "\n" );
1958 
1959  /* now print grain radius converting from cm to microns */
1960  fprintf( save.ipPnunit[ipPun], "#grn rad (mic)\tne\t" );
1961  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1962  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
1963  fprintf( save.ipPnunit[ipPun], "\n" );
1964 
1965  /* don't need to do this, ever again */
1966  lgMustPrtHeader = false;
1967  }
1968 
1969  fprintf( save.ipPnunit[ipPun], " %.5e\t%.4e",
1970  radius.depth_mid_zone ,
1971  /* electron density contributed by grains, in e/cm^3,
1972  * positive number means grain supplied free electrons */
1973  gv.TotalEden );
1974 
1975  /* average charge per grain in electrons */
1976  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1977  {
1978  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AveDustZ );
1979  }
1980  fprintf( save.ipPnunit[ipPun], "\n" );
1981  }
1982  }
1983 
1984  else if( strcmp(save.chSave[ipPun],"DUSH") == 0 )
1985  {
1986  /* grain heating */
1987  if( strcmp(chTime,"LAST") != 0 )
1988  {
1989  /* used to print header exactly one time */
1990  static bool lgMustPrtHeader = true;
1991  /* save grain charge, but do labels first if this is first zone */
1992  if( lgMustPrtHeader )
1993  {
1994  /* first print string giving grain id */
1995  fprintf( save.ipPnunit[ipPun], "#Depth" );
1996  for( size_t nd=0; nd < gv.bin.size(); ++nd )
1997  fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
1998  fprintf( save.ipPnunit[ipPun], "\n" );
1999 
2000  /* now print grain radius converting from cm to microns */
2001  fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
2002  for( size_t nd=0; nd < gv.bin.size(); ++nd )
2003  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
2004  fprintf( save.ipPnunit[ipPun], "\n" );
2005 
2006  /* don't need to do this, ever again */
2007  lgMustPrtHeader = false;
2008  }
2009  fprintf( save.ipPnunit[ipPun], " %.5e",
2010  radius.depth_mid_zone );
2011  /* grain heating */
2012  for( size_t nd=0; nd < gv.bin.size(); ++nd )
2013  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->GasHeatPhotoEl );
2014  fprintf( save.ipPnunit[ipPun], "\n" );
2015  }
2016  }
2017 
2018  else if( strcmp(save.chSave[ipPun],"DUSV") == 0 )
2019  {
2020  /* grain drift velocities */
2021  if( strcmp(chTime,"LAST") != 0 )
2022  {
2023  /* used to print header exactly one time */
2024  static bool lgMustPrtHeader = true;
2025  /* save grain velocity, but do labels first if this is first zone */
2026  if( lgMustPrtHeader )
2027  {
2028  /* first print string giving grain id */
2029  fprintf( save.ipPnunit[ipPun], "#Depth" );
2030  for( size_t nd=0; nd < gv.bin.size(); ++nd )
2031  fprintf( save.ipPnunit[ipPun], "\t%s", gv.bin[nd]->chDstLab );
2032  fprintf( save.ipPnunit[ipPun], "\n" );
2033 
2034  /* now print grain radius converting from cm to microns */
2035  fprintf( save.ipPnunit[ipPun], "#grn rad (mic)" );
2036  for( size_t nd=0; nd < gv.bin.size(); ++nd )
2037  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->AvRadius*1e4 );
2038  fprintf( save.ipPnunit[ipPun], "\n" );
2039 
2040  /* don't need to do this, ever again */
2041  lgMustPrtHeader = false;
2042  }
2043  fprintf( save.ipPnunit[ipPun], " %.5e",
2044  radius.depth_mid_zone );
2045  /* grain drift velocity in km/s */
2046  for( size_t nd=0; nd < gv.bin.size(); ++nd )
2047  fprintf( save.ipPnunit[ipPun], "\t%.3e", gv.bin[nd]->DustDftVel*1e-5 );
2048  fprintf( save.ipPnunit[ipPun], "\n" );
2049  }
2050  }
2051 
2052  /* >>chng 02 dec 30, separated scattering cross section and asymmetry factor, PvH */
2053  else if( strcmp(save.chSave[ipPun],"DUSQ") == 0 )
2054  {
2055  /* save grain Qs */
2056  if( strcmp(chTime,"LAST") == 0 )
2057  {
2058  for( j=0; j < rfield.nflux; j++ )
2059  {
2060  fprintf( save.ipPnunit[ipPun], "%11.4e",
2061  rfield.anu[j] );
2062  for( size_t nd=0; nd < gv.bin.size(); nd++ )
2063  {
2064  fprintf( save.ipPnunit[ipPun], "%9.1e%9.1e",
2065  gv.bin[nd]->dstab1[j]*4./gv.bin[nd]->IntArea,
2066  gv.bin[nd]->pure_sc1[j]*gv.bin[nd]->asym[j]*4./gv.bin[nd]->IntArea );
2067  }
2068  fprintf( save.ipPnunit[ipPun], "\n" );
2069  }
2070  }
2071  }
2072 
2073  else if( strcmp(save.chSave[ipPun],"ELEM") == 0 )
2074  {
2075  if( strcmp(chTime,"LAST") != 0 )
2076  {
2077  realnum renorm = 1.f;
2078 
2079  /* this is the index for the atomic number on the physical scale */
2080  /* >>chng 04 nov 23, use c scale throughout */
2081  nelem = (long int)save.punarg[ipPun][0];
2082  ASSERT( nelem >= ipHYDROGEN );
2083 
2084  /* don't do this if element is not turned on */
2085  if( dense.lgElmtOn[nelem] )
2086  {
2087  /* >>chng 04 nov 23, add density option, leave as cm-3
2088  * default is still norm to total of that element */
2089  if( save.punarg[ipPun][1] == 0 )
2090  renorm = dense.gas_phase[nelem];
2091 
2092  fprintf( save.ipPnunit[ipPun], " %.5e", radius.depth_mid_zone );
2093 
2094  for( j=0; j <= (nelem + 1); ++j)
2095  {
2096  fprintf( save.ipPnunit[ipPun], "\t%.2e",
2097  dense.xIonDense[nelem][j]/renorm );
2098  }
2099  if( nelem==ipHYDROGEN )
2100  {
2101  /* H2 */
2102  fprintf( save.ipPnunit[ipPun], "\t%.2e",
2103  hmi.H2_total/renorm );
2104  }
2105  /* >>chng 04 nov 23 add C and O fine structure pops */
2106  else if( nelem==ipCARBON )
2107  {
2108  fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e",
2109  colden.C1Pops[0]/renorm, colden.C1Pops[1]/renorm, colden.C1Pops[2]/renorm);
2110  fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e",
2111  colden.C2Pops[0]/renorm, colden.C2Pops[1]/renorm);
2112  fprintf( save.ipPnunit[ipPun], "\t%.2e",
2113  findspecieslocal("CO")->den/renorm );
2114  }
2115  else if( nelem==ipOXYGEN )
2116  {
2117  fprintf( save.ipPnunit[ipPun], "\t%.2e\t%.2e\t%.2e",
2118  colden.O1Pops[0]/renorm, colden.O1Pops[1]/renorm, colden.O1Pops[2]/renorm);
2119  }
2120  fprintf( save.ipPnunit[ipPun], "\n" );
2121  }
2122  }
2123  }
2124 
2125  else if( strcmp(save.chSave[ipPun],"RECA") == 0 )
2126  {
2127  /* this will create table for AGN3 then exit,
2128  * routine is in makerecom.c */
2129  ion_recombAGN( save.ipPnunit[ipPun] );
2130  cdEXIT(EXIT_FAILURE);
2131  }
2132 
2133  else if( strcmp(save.chSave[ipPun],"RECE") == 0 )
2134  {
2135  /* save recombination efficiencies,
2136  * option turned on with the "save recombination efficiencies" command
2137  * output for the save recombination coefficients command is actually
2138  * produced by a series of routines, as they generate the recombination
2139  * coefficients. these include
2140  * dielsupres, helium, hydrorecom, iibod, and makerecomb*/
2141  fprintf( save.ipPnunit[ipPun],
2142  "%12.4e %12.4e %12.4e %12.4e\n",
2143  iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].RadRecomb[ipRecRad],
2144  iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].RadRecomb[ipRecNetEsc] ,
2145  iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].RadRecomb[ipRecRad],
2146  iso_sp[ipH_LIKE][ipHYDROGEN].fb[2].RadRecomb[ipRecNetEsc]);
2147  }
2148 
2149  else if( strncmp(save.chSave[ipPun],"FEc" , 3 ) == 0 )
2150  {
2151  /* FeII continuum */
2152  if( strcmp(chTime,"LAST") == 0 )
2153  {
2154 
2155  if( save.punarg[ipPun][0] == 1 )
2156  {
2157  // row format used for grids - one time print of wavelength grid first
2158  char chTempHeader[100];
2159  long ipFeII_Cont_type;
2160  if( strcmp(save.chSave[ipPun],"FEcI") == 0 )
2161  {
2162  strcpy(chTempHeader , "#FeII inward: Wl(A)\tInt[erg cm-2 s-1]\n");
2163  ipFeII_Cont_type = 1;
2164  }
2165  else if( strcmp(save.chSave[ipPun],"FEcO") == 0 )
2166  {
2167  strcpy(chTempHeader , "#FeII outward: Wl(A)\tInt[erg cm-2 s-1]\n");
2168  ipFeII_Cont_type = 2;
2169  }
2170  else if( strcmp(save.chSave[ipPun],"FEcT") == 0 )
2171  {
2172  strcpy(chTempHeader , "#FeII total: Wl(A)\tInt[erg cm-2 s-1]\n");
2173  ipFeII_Cont_type = 3;
2174  }
2175  else
2176  TotalInsanity();
2177 
2178  if( save.lgPunHeader[ipPun] &&
2179  (!grid.lgGrid || optimize.nOptimiz==0) )
2180  {
2181  // one time print of FeII continuum header
2182  j = 0;
2183  fprintf( save.ipPnunit[ipPun], "%s%.5f",
2184  chTempHeader ,
2185  AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) ));
2186  for( j=1; j < nFeIIConBins; j++ )
2187  {
2188  fprintf( save.ipPnunit[ipPun], "\t%.5e",
2189  AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) ));
2190  }
2191  fprintf( save.ipPnunit[ipPun], "\n");
2192  save.lgPunHeader[ipPun] = false;
2193  }
2194  // give intensities
2195  fprintf( save.ipPnunit[ipPun], "%.2f",
2196  SaveFeII_cont( 0 , ipFeII_Cont_type ));
2197  for( j=1; j < nFeIIConBins; j++ )
2198  {
2199  fprintf( save.ipPnunit[ipPun], "\t%e",
2200  SaveFeII_cont( j , ipFeII_Cont_type ) );
2201  }
2202  fprintf( save.ipPnunit[ipPun], "\n");
2203  }
2204  else if( save.punarg[ipPun][0] == 2 )
2205  {
2206  // the default, four columns, wl, total, inward, outward
2207  // one time print of header
2208  if( save.lgPunHeader[ipPun] &&
2209  (!grid.lgGrid || optimize.nOptimiz==0) )
2210  {
2211  fprintf( save.ipPnunit[ipPun], "FeiI wl, total, inward, outward\n");
2212  save.lgPunHeader[ipPun] = false;
2213  }
2214  for( j=0; j < nFeIIConBins; j++ )
2215  {
2216  fprintf( save.ipPnunit[ipPun], "%.5f",
2217  AnuUnit((realnum)RYDLAM/((FeII_Cont[j][0]+FeII_Cont[j+1][0])/2.f) ));
2218  fprintf( save.ipPnunit[ipPun], "\t%e",
2219  SaveFeII_cont( j , 3 ) );// total
2220  fprintf( save.ipPnunit[ipPun], "\t%e",
2221  SaveFeII_cont( j , 1 ) );// inward
2222  fprintf( save.ipPnunit[ipPun], "\t%e",
2223  SaveFeII_cont( j , 2 ) );// outward
2224  fprintf( save.ipPnunit[ipPun], "\n");
2225  }
2226  }
2227  }
2228  }
2229 
2230  /* save column densities */
2231  else if( strcmp(save.chSave[ipPun],"FENl") == 0 )
2232  {
2233  if( strcmp(chTime,"LAST") == 0 )
2234  {
2235  /* save FeII level energies and stat weights, followed by column density */
2236  FeIIPunchColden( save.ipPnunit[ipPun] );
2237  }
2238  }
2239 
2240  else if( strcmp(save.chSave[ipPun],"FE2l") == 0 )
2241  {
2242  if( strcmp(chTime,"LAST") == 0 )
2243  {
2244  /* save FeII level energies and stat weights */
2245  FeIIPunchLevels( save.ipPnunit[ipPun] );
2246  }
2247  }
2248 
2249  else if( strcmp(save.chSave[ipPun],"FEli") == 0 )
2250  {
2251  if( strcmp(chTime,"LAST") == 0 )
2252  {
2253  /* save line intensities, routine is in atom_feii.c */
2254  FeIISaveLines( save.ipPnunit[ipPun] );
2255  }
2256  }
2257 
2258  else if( strcmp(save.chSave[ipPun],"FE2o") == 0 )
2259  {
2260  if( strcmp(chTime,"LAST") == 0 )
2261  {
2262  /* save line optical depths, routine is in atom_feii.c */
2263  FeIIPunchOpticalDepth( save.ipPnunit[ipPun] );
2264  }
2265  }
2266 
2267  else if( strcmp(save.chSave[ipPun],"FRED") == 0 )
2268  {
2269  /* set with save Fred command, this punches some stuff from
2270  * Fred Hamann's dynamics project */
2271  if( strcmp(chTime,"LAST") != 0 )
2272  {
2273  /* Fred's list */
2274  fprintf( save.ipPnunit[ipPun], "%.5e\t%.5e\t%.3e\t%.3e\t%.3e"
2275  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
2276  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
2277  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
2278  "\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e"
2279  "\t%.3e\t%.3e\n",
2280  radius.Radius, radius.depth ,wind.windv/1e5,
2281  wind.dvdr,
2282  dense.gas_phase[ipHYDROGEN], dense.eden , phycon.te,
2283  wind.AccelLine , wind.AccelCont ,
2284  wind.fmul ,
2285  // acceleration in this zone due to electron scattering,
2286  // if incident SED was not attenuated
2287  pressure.pinzon_PresIntegElecThin/dense.xMassDensity/radius.drad_x_fillfac ,
2288  mean.xIonMean[0][ipHYDROGEN][0][0] , mean.xIonMean[0][ipHYDROGEN][1][0] ,
2289  mean.xIonMean[0][ipHELIUM][0][0] , mean.xIonMean[0][ipHELIUM][1][0] ,
2290  mean.xIonMean[0][ipHELIUM][2][0] ,
2291  mean.xIonMean[0][ipCARBON][1][0] , mean.xIonMean[0][ipCARBON][2][0] ,
2292  mean.xIonMean[0][ipCARBON][3][0] ,
2293  mean.xIonMean[0][ipOXYGEN][0][0] , mean.xIonMean[0][ipOXYGEN][1][0] ,
2294  mean.xIonMean[0][ipOXYGEN][2][0] , mean.xIonMean[0][ipOXYGEN][3][0] ,
2295  mean.xIonMean[0][ipOXYGEN][4][0] , mean.xIonMean[0][ipOXYGEN][5][0] ,
2296  mean.xIonMean[0][ipOXYGEN][6][0] , mean.xIonMean[0][ipOXYGEN][7][0] ,
2297  dense.xIonDense[ipHYDROGEN][0] , dense.xIonDense[ipHYDROGEN][1] ,
2298  dense.xIonDense[ipHELIUM][0] , dense.xIonDense[ipHELIUM][1] ,
2299  dense.xIonDense[ipHELIUM][2] ,
2300  dense.xIonDense[ipCARBON][1] , dense.xIonDense[ipCARBON][2] ,
2301  dense.xIonDense[ipCARBON][3] ,
2302  dense.xIonDense[ipOXYGEN][0] , dense.xIonDense[ipOXYGEN][1] ,
2303  dense.xIonDense[ipOXYGEN][2] , dense.xIonDense[ipOXYGEN][3] ,
2304  dense.xIonDense[ipOXYGEN][4] , dense.xIonDense[ipOXYGEN][5] ,
2305  dense.xIonDense[ipOXYGEN][6] , dense.xIonDense[ipOXYGEN][7] ,
2306  mean.xIonMean[0][ipMAGNESIUM][1][0] , dense.xIonDense[ipMAGNESIUM][1],
2307  TauLines[ipT1032].Emis().TauIn() ,
2308  TauLines[ipT1032].Emis().TauCon() );
2309  }
2310  }
2311 
2312  else if( strcmp(save.chSave[ipPun],"FE2d") == 0 )
2313  {
2314  /* save some departure coefficients for large FeII atom */
2315  if( strcmp(chTime,"LAST") != 0 )
2316  FeIIPunDepart(save.ipPnunit[ipPun],false);
2317  }
2318 
2319  else if( strcmp(save.chSave[ipPun],"FE2D") == 0 )
2320  {
2321  /* save all departure coefficients for large FeII atom */
2322  if( strcmp(chTime,"LAST") != 0 )
2323  FeIIPunDepart(save.ipPnunit[ipPun],true);
2324  }
2325 
2326  else if( strcmp(save.chSave[ipPun],"FE2p") == 0 )
2327  {
2328  bool lgFlag = false;
2329  if( save.punarg[ipPun][2] )
2330  lgFlag = true;
2331  /* save small subset of level populations for large FeII atom */
2332  if( strcmp(chTime,"LAST") != 0 )
2333  FeIIPunPop(save.ipPnunit[ipPun],false,0,0,
2334  lgFlag);
2335  }
2336 
2337  else if( strcmp(save.chSave[ipPun],"FE2P") == 0 )
2338  {
2339  bool lgFlag = false;
2340  if( save.punarg[ipPun][2] )
2341  lgFlag = true;
2342  /* save range of level populations for large FeII atom */
2343  if( strcmp(chTime,"LAST") != 0 )
2344  FeIIPunPop(save.ipPnunit[ipPun],
2345  true,
2346  (long int)save.punarg[ipPun][0],
2347  (long int)save.punarg[ipPun][1],
2348  lgFlag);
2349  }
2350 
2351  /* save spectra in fits format */
2352  else if( strcmp(save.chSave[ipPun],"FITS") == 0 )
2353  {
2354  if( strcmp(chTime,"LAST") == 0 )
2355  {
2356  saveFITSfile( save.ipPnunit[ipPun], NUM_OUTPUT_TYPES );
2357  }
2358  }
2359  /* save gammas (but without element) */
2360  else if( strcmp(save.chSave[ipPun],"GAMt") == 0 )
2361  {
2362  if( strcmp(chTime,"LAST") != 0 )
2363  {
2364  long ns;
2365  /* save photoionization rates, with the PUNCH GAMMAS command */
2366  for( nelem=0; nelem < LIMELM; nelem++ )
2367  {
2368  for( ion=0; ion <= nelem; ion++ )
2369  {
2370  for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
2371  {
2372  fprintf( save.ipPnunit[ipPun], "%3ld%3ld%3ld%10.2e%10.2e%10.2e",
2373  nelem+1, ion+1, ns+1,
2374  ionbal.PhotoRate_Shell[nelem][ion][ns][0],
2375  ionbal.PhotoRate_Shell[nelem][ion][ns][1] ,
2376  ionbal.PhotoRate_Shell[nelem][ion][ns][2] );
2377 
2378  for( j=0; j < t_yield::Inst().nelec_eject(nelem,ion,ns); j++ )
2379  {
2380  fprintf( save.ipPnunit[ipPun], "%5.2f",
2381  t_yield::Inst().elec_eject_frac(nelem,ion,ns,j) );
2382  }
2383  fprintf( save.ipPnunit[ipPun], "\n" );
2384  }
2385  }
2386  }
2387  }
2388  }
2389 
2390  /* save gammas element, ion */
2391  else if( strcmp(save.chSave[ipPun],"GAMe") == 0 )
2392  {
2393  if( strcmp(chTime,"LAST") != 0 )
2394  {
2395  int ns;
2396  nelem = (long)save.punarg[ipPun][0];
2397  ion = (long)save.punarg[ipPun][1];
2398  /* valence shell */
2399  ns = Heavy.nsShells[nelem][ion]-1;
2400  /* show what some of the ionization sources are */
2401  GammaPrt(
2402  opac.ipElement[nelem][ion][ns][0] ,
2403  opac.ipElement[nelem][ion][ns][1] ,
2404  opac.ipElement[nelem][ion][ns][2] ,
2405  save.ipPnunit[ipPun],
2406  ionbal.PhotoRate_Shell[nelem][ion][ns][0] ,
2407  ionbal.PhotoRate_Shell[nelem][ion][ns][0]*0.1 );
2408  }
2409  }
2410 
2411  else if( strcmp(save.chSave[ipPun],"GAUN") == 0 )
2412  {
2413  /* save gaunt factors */
2414  if( strcmp(chTime,"LAST") != 0 )
2415  SaveGaunts(save.ipPnunit[ipPun]);
2416  }
2417 
2418  // generating the SAVE GRID output has been moved to cdPrepareExit()
2419  // to make sure that the output always records any type of failure
2420  }
2421  else
2422  {
2423  //no hit this branch, key should be in next
2424  lgNoHitFirstBranch = true;
2425  }
2426 
2427  // hack needed for code to compile with Visual Studio
2428  // keep this identical to the if-statement further up!!
2429  if( iterations.lgLastIt || !save.lgPunLstIter[ipPun] ||
2430  ( lgAbort && strcmp(chTime,"LAST") == 0 ) ||
2431  ( dynamics.lgTimeDependentStatic && dynamics.lgStatic_completed ) )
2432  {
2433  if( strcmp(save.chSave[ipPun],"HISp") == 0 )
2434  {
2435  /* save pressure history of current zone */
2436  if( strcmp(chTime,"LAST") != 0 )
2437  {
2438  /* note if pressure convergence failure occurred in history that follows */
2439  if( !conv.lgConvPres )
2440  {
2441  fprintf( save.ipPnunit[ipPun],
2442  "#PROBLEM Pressure not converged iter %li zone %li density-pressure follows:\n",
2443  iteration , nzone );
2444  }
2445  /* note if temperature convergence failure occurred in history that follows */
2446  if( !conv.lgConvTemp )
2447  {
2448  fprintf( save.ipPnunit[ipPun],
2449  "#PROBLEM Temperature not converged iter %li zone %li density-pressure follows:\n",
2450  iteration , nzone );
2451  }
2452  for( unsigned long k=0; k < conv.hist_pres_density.size(); ++k )
2453  {
2454  /* save history of density - pressure, with correct pressure */
2455  fprintf( save.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
2456  iteration,
2457  nzone,
2458  conv.hist_pres_density[k],
2459  conv.hist_pres_current[k],
2460  conv.hist_pres_error[k]);
2461  }
2462  }
2463  }
2464 
2465  else if( strcmp(save.chSave[ipPun],"HISt") == 0 )
2466  {
2467  /* save temperature history of current zone */
2468  if( strcmp(chTime,"LAST") != 0 )
2469  {
2470  /* note if pressure convergence failure occurred in history that follows */
2471  if( !conv.lgConvPres )
2472  {
2473  fprintf( save.ipPnunit[ipPun],
2474  "#PROBLEM Pressure not converged iter %li zone %li temp heat cool follows:\n",
2475  iteration , nzone );
2476  }
2477  /* note if temperature convergence failure occurred in history that follows */
2478  if( !conv.lgConvTemp )
2479  {
2480  fprintf( save.ipPnunit[ipPun],
2481  "#PROBLEM Temperature not converged iter %li zone %li temp heat cool follows:\n",
2482  iteration , nzone );
2483  }
2484  for( unsigned long k=0; k < conv.hist_temp_temp.size(); ++k )
2485  {
2486  /* save history of density - pressure, with correct pressure */
2487  fprintf( save.ipPnunit[ipPun] , "%2li %4li\t%.5e\t%.5e\t%.5e\n",
2488  iteration,
2489  nzone,
2490  conv.hist_temp_temp[k],
2491  conv.hist_temp_heat[k],
2492  conv.hist_temp_cool[k]);
2493  }
2494  }
2495  }
2496 
2497  else if( strncmp(save.chSave[ipPun],"H2",2) == 0 )
2498  {
2499  /* all save info on large H2 molecule include H2 PDR pdr */
2500  save.whichDiatomToPrint[ipPun]->H2_PunchDo( save.ipPnunit[ipPun] , save.chSave[ipPun] , chTime, ipPun );
2501  }
2502 
2503  else if( strcmp(save.chSave[ipPun],"HEAT") == 0 )
2504  {
2505  /* save heating, routine in file of same name */
2506  if( strcmp(chTime,"LAST") != 0 )
2507  SaveHeat(save.ipPnunit[ipPun]);
2508  }
2509 
2510  else if( strncmp(save.chSave[ipPun],"HE",2) == 0 )
2511  {
2512  /* various save helium commands */
2513  /* save helium line wavelengths */
2514  if( strcmp(save.chSave[ipPun] , "HELW") == 0 )
2515  {
2516  if( strcmp(chTime,"LAST") == 0 )
2517  {
2518  /* save helium & he-like wavelengths, first header */
2519  fprintf( save.ipPnunit[ipPun],
2520  "Z\tElem\t2 1P->1 1S\t2 3P1->1 1S\t2 3P2->1 1S"
2521  "\t2 3S->1 1S\t2 3P2->2 3S\t2 3P1->2 3S\t2 3P0->2 3S" );
2522  fprintf( save.ipPnunit[ipPun], "\n" );
2523  for( nelem=ipHELIUM; nelem<LIMELM; ++nelem )
2524  {
2525  /* print element name, nuclear charge */
2526  fprintf( save.ipPnunit[ipPun], "%li\t%s",
2527  nelem+1 , elementnames.chElementSym[nelem] );
2528  /*prt_wl print floating wavelength in Angstroms, in output format */
2529  fprintf( save.ipPnunit[ipPun], "\t" );
2530  prt_wl( save.ipPnunit[ipPun] ,
2531  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p1P,ipHe1s1S).WLAng() );
2532  fprintf( save.ipPnunit[ipPun], "\t" );
2533  prt_wl( save.ipPnunit[ipPun] ,
2534  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P1,ipHe1s1S).WLAng() );
2535  fprintf( save.ipPnunit[ipPun], "\t" );
2536  prt_wl( save.ipPnunit[ipPun] ,
2537  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P2,ipHe1s1S).WLAng() );
2538  fprintf( save.ipPnunit[ipPun], "\t" );
2539  prt_wl( save.ipPnunit[ipPun] ,
2540  iso_sp[ipHE_LIKE][nelem].trans(ipHe2s3S,ipHe1s1S).WLAng() );
2541  fprintf( save.ipPnunit[ipPun], "\t" );
2542  prt_wl( save.ipPnunit[ipPun] ,
2543  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P2,ipHe2s3S).WLAng() );
2544  fprintf( save.ipPnunit[ipPun], "\t" );
2545  prt_wl( save.ipPnunit[ipPun] ,
2546  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P1,ipHe2s3S).WLAng() );
2547  fprintf( save.ipPnunit[ipPun], "\t" );
2548  prt_wl( save.ipPnunit[ipPun] ,
2549  iso_sp[ipHE_LIKE][nelem].trans(ipHe2p3P0,ipHe2s3S).WLAng() );
2550  fprintf( save.ipPnunit[ipPun], "\n");
2551  }
2552  }
2553  }
2554  else
2555  TotalInsanity();
2556  }
2557 
2558  /* save hummer, results needed for Lya transport, to feed into David's routine */
2559  else if( strcmp(save.chSave[ipPun],"HUMM") == 0 )
2560  {
2561  eps = iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Aul()/
2562  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Coll().ColUL( colliders );
2563  fprintf( save.ipPnunit[ipPun],
2564  " %.5e %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n",
2565  radius.depth_mid_zone,
2566  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn(),
2567  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop(),
2568  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop(),
2569  phycon.te, iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().damp(), eps );
2570  }
2571 
2572  else if( strncmp( save.chSave[ipPun] , "HYD", 3 ) == 0 )
2573  {
2574  /* various save hydrogen commands */
2575  if( strcmp(save.chSave[ipPun],"HYDc") == 0 )
2576  {
2577  if( strcmp(chTime,"LAST") != 0 )
2578  {
2579  /* save hydrogen physical conditions */
2580  fprintf( save.ipPnunit[ipPun],
2581  " %.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
2582  radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN], dense.eden,
2583  dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
2584  dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN],
2585  hmi.H2_total/dense.gas_phase[ipHYDROGEN],
2586  findspecieslocal("H2+")->den/dense.gas_phase[ipHYDROGEN],
2587  findspecieslocal("H3+")->den/dense.gas_phase[ipHYDROGEN],
2588  findspecieslocal("H-")->den/dense.gas_phase[ipHYDROGEN] );
2589  }
2590  }
2591 
2592  else if( strcmp(save.chSave[ipPun],"HYDi") == 0 )
2593  {
2594  if( strcmp(chTime,"LAST") != 0 )
2595  {
2596  /* save hydrogen ionization
2597  * this will be total decays to ground */
2598  RateInter = 0.;
2599  stage = iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ColIoniz*dense.eden*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
2600  fref = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
2601  fout = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop();
2602  /* 06 aug 28, from numLevels_max to _local. */
2603  for( ion=ipH2s; ion < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local; ion++ )
2604  {
2605  /* this is total decays to ground */
2606  RateInter +=
2607  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Aul()*
2608  (iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pesc() +
2609  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pelec_esc() +
2610  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ion,ipH1s).Emis().Pdest());
2611  /* total photo from all levels */
2612  fref += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ion].gamnc*iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
2613  /* total col ion from all levels */
2614  stage += iso_sp[ipH_LIKE][ipHYDROGEN].fb[ion].ColIoniz*dense.eden*
2615  iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
2616  fout += iso_sp[ipH_LIKE][ipHYDROGEN].st[ion].Pop();
2617  }
2618 
2619  /* make these relative to parent ion */
2620  stage /= dense.xIonDense[ipHYDROGEN][1];
2621  fref /= dense.xIonDense[ipHYDROGEN][1];
2622  fout /= dense.xIonDense[ipHYDROGEN][1];
2623 
2624  fprintf( save.ipPnunit[ipPun], "hion\t%4ld\t%.2e\t%.2e\t%.2e",
2625  nzone,
2626  /* photo and collision ion rates have units s-1 */
2627  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc,
2628  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ColIoniz* dense.EdenHCorr,
2629  ionbal.RateRecomTot[ipHYDROGEN][0] );
2630 
2631  fprintf( save.ipPnunit[ipPun], "\t%.2e",
2632  iso_sp[ipH_LIKE][ipHYDROGEN].RadRec_caseB );
2633 
2634  fprintf( save.ipPnunit[ipPun],
2635  "\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
2636  dense.xIonDense[ipHYDROGEN][1]/dense.xIonDense[ipHYDROGEN][0],
2637  iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc/(ionbal.RateRecomTot[ipHYDROGEN][0]),
2638  iso_sp[ipH_LIKE][ipHYDROGEN].fb[1].RadRecomb[ipRecEsc],
2639  RateInter,
2640  fref/MAX2(1e-37,fout),
2641  stage/MAX2(1e-37,fout),
2642  /* simple H+ */
2643  safe_div( iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*dense.xIonDense[ipHYDROGEN][0], dense.eden*dense.xIonDense[ipHYDROGEN][1] ),
2644  secondaries.csupra[ipHYDROGEN][0]);
2645 
2646  GammaPrt(iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon,rfield.nflux,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipOpac,
2647  save.ipPnunit[ipPun],iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc,iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].gamnc*
2648  0.05);
2649  }
2650  }
2651 
2652  else if( strcmp(save.chSave[ipPun],"HYDp") == 0 )
2653  {
2654  if( strcmp(chTime,"LAST") != 0 )
2655  {
2656  /* save hydrogen populations
2657  * first give total atom and ion density [cm-3]*/
2658  fprintf( save.ipPnunit[ipPun], "%.5e\t%.2e\t%.2e",
2659  radius.depth_mid_zone,
2660  dense.xIonDense[ipHYDROGEN][0],
2661  dense.xIonDense[ipHYDROGEN][1] );
2662 
2663  /* next give state-specific densities [cm-3] */
2664  for( j=ipH1s; j < iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local-1; j++ )
2665  {
2666  fprintf( save.ipPnunit[ipPun], "\t%.2e",
2667  iso_sp[ipH_LIKE][ipHYDROGEN].st[j].Pop() );
2668  }
2669  fprintf( save.ipPnunit[ipPun], "\n" );
2670  }
2671  }
2672 
2673  else if( strcmp(save.chSave[ipPun],"HYDl") == 0 )
2674  {
2675  if( strcmp(chTime,"LAST") == 0 )
2676  {
2677  /* save hydrogen line
2678  * gives intensities and optical depths */
2679  for( ipHi=1; ipHi<iso_sp[ipH_LIKE][ipHYDROGEN].numLevels_local -
2680  iso_sp[ipH_LIKE][ipHYDROGEN].nCollapsed_local; ++ipHi )
2681  {
2682  for( ipLo=0; ipLo<ipHi; ++ipLo )
2683  {
2684  if( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).ipCont() < 0 )
2685  continue;
2686  fprintf(save.ipPnunit[ipPun], "%li\t%li\t%li\t%li\t%.4e\t%.2e\n",
2687  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipHi].n(),
2688  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipHi].l(),
2689  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipLo].n(),
2690  iso_sp[ipH_LIKE][ipHYDROGEN].st[ipLo].l(),
2691  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).EnergyRyd(),
2692  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipHi,ipLo).Emis().TauIn() );
2693  }
2694  }
2695  }
2696  }
2697 
2698  /* save hydrogen Lya - some details about Lya */
2699  else if( strcmp(save.chSave[ipPun],"HYDL") == 0 )
2700  {
2701  if( strcmp(chTime,"LAST") != 0 )
2702  {
2703  /* the population ratio for Lya */
2704  double popul = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH2p].Pop()/SDIV(iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop());
2705  /* the excitation temperature of Lya */
2706  texc = TexcLine( iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s) );
2707  fprintf( save.ipPnunit[ipPun],
2708  "%.5e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
2709  radius.depth_mid_zone,
2710  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauIn(),
2711  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().TauTot(),
2712  popul,
2713  texc,
2714  phycon.te,
2715  texc/phycon.te ,
2716  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pesc(),
2717  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().Pdest(),
2718  iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).Emis().pump(),
2719  opac.opacity_abs[iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont()-1],
2720  opac.albedo[iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,ipH1s).ipCont()-1] );
2721  }
2722  }
2723 
2724  else if( strcmp(save.chSave[ipPun],"HYDr") == 0 )
2725  {
2726  /* save hydrogen recc - recombination cooling for AGN3 */
2727  TempChange(2500.f, false);
2728  while( phycon.te <= 20000. )
2729  {
2730  double r1;
2731  double ThinCoolingCaseB;
2732 
2733  r1 = HydroRecCool(1,0);
2734  ThinCoolingCaseB = pow(10.,((-25.859117 +
2735  0.16229407*phycon.telogn[0] +
2736  0.34912863*phycon.telogn[1] -
2737  0.10615964*phycon.telogn[2])/(1. +
2738  0.050866793*phycon.telogn[0] -
2739  0.014118924*phycon.telogn[1] +
2740  0.0044980897*phycon.telogn[2] +
2741  6.0969594e-5*phycon.telogn[3])))/phycon.te;
2742 
2743  fprintf( save.ipPnunit[ipPun], " %10.2e\t",
2744  phycon.te);
2745  fprintf( save.ipPnunit[ipPun], " %10.2e\t",
2746  (r1+ThinCoolingCaseB)/(BOLTZMANN*phycon.te) );
2747 
2748  fprintf( save.ipPnunit[ipPun], " %10.2e\t",
2749  r1/(BOLTZMANN*phycon.te));
2750 
2751  fprintf( save.ipPnunit[ipPun], " %10.2e\n",
2752  ThinCoolingCaseB/(BOLTZMANN*phycon.te));
2753 
2754  TempChange(phycon.te *2.f , false);
2755  }
2756  /* must exit since we have disturbed the solution */
2757  fprintf(ioQQQ , "save agn now exits since solution is disturbed.\n");
2758  cdEXIT( EXIT_SUCCESS );
2759  }
2760  else
2761  TotalInsanity();
2762  }
2763 
2764  else if( strcmp(save.chSave[ipPun],"IONI") == 0 )
2765  {
2766  if( strcmp(chTime,"LAST") == 0 )
2767  {
2768  /* save mean ionization distribution */
2769  PrtMeanIon( 'i', false , save.ipPnunit[ipPun] );
2770  }
2771  }
2772 
2773  /* save ionization rates */
2774  else if( strcmp(save.chSave[ipPun],"IONR") == 0 )
2775  {
2776  if( strcmp(chTime,"LAST") != 0 )
2777  {
2778  /* this is element number */
2779  nelem = (long)save.punarg[ipPun][0];
2780  fprintf( save.ipPnunit[ipPun],
2781  "%.5e\t%.4e\t%.4e",
2782  radius.depth_mid_zone,
2783  dense.eden ,
2784  dynamics.Rate);
2785  /* >>chng 04 oct 15, from nelem+2 to nelem+1 - array over read -
2786  * caught by PnH */
2787  for( ion=0; ion<nelem+1; ++ion )
2788  {
2789  fprintf( save.ipPnunit[ipPun],
2790  "\t%.4e\t%.4e\t%.4e\t%.4e",
2791  dense.xIonDense[nelem][ion] ,
2792  ionbal.RateIonizTot(nelem,ion) ,
2793  ionbal.RateRecomTot[nelem][ion] ,
2794  dynamics.Source[nelem][ion] );
2795  }
2796  fprintf( save.ipPnunit[ipPun], "\n");
2797  }
2798  }
2799 
2800  else if( strcmp(save.chSave[ipPun]," IP ") == 0 )
2801  {
2802  if( strcmp(chTime,"LAST") == 0 )
2803  {
2804  /* save valence shell ip's */
2805  for( nelem=0; nelem < LIMELM; nelem++ )
2806  {
2807  int ion_big;
2808  double energy;
2809 
2810  /* this is the largest number of ion stages per line */
2811  const int NELEM_LINE = 10;
2812  /* this loop in case all ions do not fit across page */
2813  for( ion_big=0; ion_big<=nelem; ion_big += NELEM_LINE )
2814  {
2815  int ion_limit = MIN2(ion_big+NELEM_LINE-1,nelem);
2816 
2817  /* new line then element name */
2818  fprintf( save.ipPnunit[ipPun],
2819  "\n%2.2s", elementnames.chElementSym[nelem]);
2820 
2821  /* print ion stages across line */
2822  for( ion=ion_big; ion <= ion_limit; ++ion )
2823  {
2824  fprintf( save.ipPnunit[ipPun], "\t%4ld", ion+1 );
2825  }
2826  fprintf( save.ipPnunit[ipPun], "\n" );
2827 
2828  /* this loop is over all shells */
2829  ASSERT( ion_limit < LIMELM );
2830  /* upper limit is number of shells in atom */
2831  for( ips=0; ips < Heavy.nsShells[nelem][ion_big]; ips++ )
2832  {
2833 
2834  /* print shell label */
2835  fprintf( save.ipPnunit[ipPun], "%2.2s", Heavy.chShell[ips]);
2836 
2837  /* loop over possible ions */
2838  for( ion=ion_big; ion<=ion_limit; ++ion )
2839  {
2840 
2841  /* does this subshell exist for this ion? break if it does not*/
2842  /*if( Heavy.nsShells[nelem][ion]<Heavy.nsShells[nelem][0] )*/
2843  if( ips >= Heavy.nsShells[nelem][ion] )
2844  break;
2845 
2846  /* array elements are shell, numb of electrons, element, 0 */
2847  energy = t_ADfA::Inst().ph1(ips,nelem-ion,nelem,0);
2848 
2849  /* now print threshold with correct format */
2850  if( energy < 10. )
2851  {
2852  fprintf( save.ipPnunit[ipPun], "\t%6.3f", energy );
2853  }
2854  else if( energy < 100. )
2855  {
2856  fprintf( save.ipPnunit[ipPun], "\t%6.2f", energy );
2857  }
2858  else if( energy < 1000. )
2859  {
2860  fprintf( save.ipPnunit[ipPun], "\t%6.1f", energy );
2861  }
2862  else
2863  {
2864  fprintf( save.ipPnunit[ipPun], "\t%6ld", (long)(energy) );
2865  }
2866  }
2867 
2868  /* put cs at end of long line */
2869  fprintf( save.ipPnunit[ipPun], "\n" );
2870  }
2871  }
2872  }
2873  }
2874  }
2875 
2876  else if( strcmp(save.chSave[ipPun],"LINC") == 0 )
2877  {
2878  /* save line cumulative */
2879  if( strcmp(chTime,"LAST") != 0 )
2880  {
2881  save_line(save.ipPnunit[ipPun],"PUNC",
2882  save.lgEmergent[ipPun]);
2883  }
2884  }
2885 
2886  else if( strcmp(save.chSave[ipPun],"LIND") == 0 )
2887  {
2888  /* save line data, then stop */
2889  SaveLineData(save.ipPnunit[ipPun]);
2890  }
2891 
2892  else if( strcmp(save.chSave[ipPun],"LINL") == 0 )
2893  {
2894  /* save line labels */
2895  bool lgPrintAll=false;
2896  /* LONG keyword on save line labels command sets this to 1 */
2897  if( save.punarg[ipPun][0]>0. )
2898  lgPrintAll = true;
2899  prt_LineLabels(save.ipPnunit[ipPun] , lgPrintAll );
2900  }
2901 
2902  else if( strcmp(save.chSave[ipPun],"LINO") == 0 )
2903  {
2904  if( strcmp(chTime,"LAST") == 0 )
2905  {
2906  /* save line optical depths, routine is below, file static */
2907  SaveLineStuff(save.ipPnunit[ipPun],"optical" , save.punarg[ipPun][0]);
2908  }
2909  }
2910 
2911  else if( strcmp(save.chSave[ipPun],"LINP") == 0 )
2912  {
2913  if( strcmp(chTime,"LAST") != 0 )
2914  {
2915  static bool lgFirst=true;
2916  /* save line populations, need to do this twice if very first
2917  * call since first call to SaveLineStuff generates atomic parameters
2918  * rather than level pops, routine is below, file static */
2919  SaveLineStuff(save.ipPnunit[ipPun],"populat" , save.punarg[ipPun][0]);
2920  if( lgFirst )
2921  {
2922  lgFirst = false;
2923  SaveLineStuff(save.ipPnunit[ipPun],"populat" , save.punarg[ipPun][0]);
2924  }
2925  }
2926  }
2927 
2928  else if( strcmp(save.chSave[ipPun],"LINS") == 0 )
2929  {
2930  /* save line emissivity */
2931  if( strcmp(chTime,"LAST") != 0 )
2932  {
2933  save_line(save.ipPnunit[ipPun],"PUNS",
2934  save.lgEmergent[ipPun]);
2935  }
2936  }
2937 
2938  else if( strcmp(save.chSave[ipPun],"LINR") == 0 )
2939  {
2940  /* save line RT */
2941  if( strcmp(chTime,"LAST") != 0 )
2942  Save_Line_RT( save.ipPnunit[ipPun]);
2943  }
2944 
2945  else if( strcmp(save.chSave[ipPun],"LINA") == 0 )
2946  {
2947  /* save line array */
2948  if( strcmp(chTime,"LAST") == 0 )
2949  {
2950  /* save out all lines with energies */
2951  for( j=0; j < LineSave.nsum; j++ )
2952  {
2953  if( LineSv[j].wavelength > 0. &&
2954  LineSv[j].SumLine[0] > 0. )
2955  {
2956  /* line energy, in units set with units option */
2957  fprintf( save.ipPnunit[ipPun], "%12.5e",
2958  AnuUnit((realnum)RYDLAM/LineSv[j].wavelength) );
2959  /* line label */
2960  fprintf( save.ipPnunit[ipPun], "\t%4.4s ",
2961  LineSv[j].chALab );
2962  /* wavelength */
2963  prt_wl( save.ipPnunit[ipPun], LineSv[j].wavelength );
2964  /* intrinsic intensity */
2965  fprintf( save.ipPnunit[ipPun], "\t%8.3f",
2966  log10(SDIV(LineSv[j].SumLine[0]) ) + radius.Conv2PrtInten );
2967  /* emergent line intensity, r recombination */
2968  fprintf( save.ipPnunit[ipPun], "\t%8.3f",
2969  log10(SDIV(LineSv[j].SumLine[1]) ) + radius.Conv2PrtInten );
2970  /* type of line, i for info, etc */
2971  fprintf( save.ipPnunit[ipPun], " \t%c\n",
2972  LineSv[j].chSumTyp);
2973  }
2974  }
2975  }
2976  }
2977 
2978  else if( strcmp(save.chSave[ipPun],"LINI") == 0 )
2979  {
2980  if( strcmp(chTime,"LAST") == 0 &&
2981  (nzone/save.LinEvery)*save.LinEvery != nzone )
2982  {
2983  /* this is the last zone
2984  * save line intensities - but do not do last zone twice */
2985  SaveLineIntensity(save.ipPnunit[ipPun] , ipPun , save.punarg[ipPun][0] );
2986  }
2987  else if( strcmp(chTime,"LAST") != 0 )
2988  {
2989  /* following so we only save first zone if LinEvery reset */
2990  if( (save.lgLinEvery && nzone == 1) ||
2991  (nzone/save.LinEvery)*save.LinEvery == nzone )
2992  {
2993  /* this is middle of calculation
2994  * save line intensities */
2995  SaveLineIntensity(save.ipPnunit[ipPun] , ipPun , save.punarg[ipPun][0]);
2996  }
2997  }
2998  }
2999 
3000  else if( strcmp( save.chSave[ipPun],"LEIL") == 0)
3001  {
3002  /* some line intensities for the Leiden PDR,
3003  * but only do this when calculation is complete */
3004  if( strcmp(chTime,"LAST") == 0 )
3005  {
3006  double absval , rel;
3007  long int n;
3008  /* the lines we will find,
3009  * for a sample list of PDR lines look at LineList_PDR_H2.dat
3010  * in the cloudy data dir */
3011  /* the number of H2 lines */
3012  const int NLINE_H2 = 31;
3013  /* the number of lines which are not H2 */
3014  const int NLINE_NOTH_H2 = 5;
3015  /* the labels and wavelengths for the lines that are not H2 */
3016  char chLabel[NLINE_NOTH_H2][5]=
3017  {"C 2", "O 1", "O 1","C 1", "C 1" };
3018  double Wl[NLINE_NOTH_H2]=
3019  {157.6 , 63.17 , 145.5 ,609.2 , 369.7 , };
3020  /* these are wavelengths in microns, conv to Angstroms before call */
3021  /* >>chng 05 sep 06, many of following wavelengths updated to agree
3022  * with output - apparently not updated when energies changed */
3023  double Wl_H2[NLINE_H2]=
3024  {2.121 ,
3025  28.213, 17.03 , 12.28 , 9.662 , 8.024 , 6.907 , 6.107 , 5.510 , 5.051 , 4.693 ,
3026  4.408 , 4.180 , 3.996 , 3.845 , 3.723 , 3.625 , 3.546 , 3.485 , 3.437 , 3.403 ,
3027  3.380 , 3.368 , 3.365 , 3.371 , 3.387 , 3.410 , 3.441 , 3.485 , 3.542 , 3.603};
3028  /* print a header for the lines */
3029  for( n=0; n<NLINE_NOTH_H2; ++n )
3030  {
3031  fprintf(save.ipPnunit[ipPun], "%s\t%.2f",chLabel[n] , Wl[n]);
3032  /* get the line, non positive return says didn't find it */
3033  /* arguments are 4-char label, wavelength, return log total intensity, linear rel inten */
3034  if( cdLine( chLabel[n] , (realnum)(Wl[n]*1e4) , &absval , &rel ) <= 0 )
3035  {
3036  fprintf(save.ipPnunit[ipPun], " did not find\n");
3037  }
3038  else
3039  {
3040  fprintf(save.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
3041  }
3042  }
3043  fprintf(save.ipPnunit[ipPun], "\n\n\n");
3044 
3045  /* only print the H2 lines if the big molecule is turned on */
3046  if( h2.lgEnabled )
3047  {
3048  fprintf(save.ipPnunit[ipPun],
3049  "Here are some of the H2 Intensities, The first one is the\n"
3050  "1-0 S(0) line and the following ones are the 0-0 S(X)\n"
3051  "lines where X goes from 0 to 29\n\n");
3052  for( n=0; n<NLINE_H2; ++n )
3053  {
3054  fprintf(save.ipPnunit[ipPun], "%s\t%.3f","H2 " , Wl_H2[n]);
3055  /* get the line, non positive return says didn't find it */
3056  if( cdLine( "H2 " , (realnum)(Wl_H2[n]*1e4) , &absval , &rel ) <= 0 )
3057  {
3058  fprintf(save.ipPnunit[ipPun], " did not find\n");
3059  }
3060  else
3061  {
3062  fprintf(save.ipPnunit[ipPun], "\t%.3e\t%.3e\n",pow(10.,rel),absval);
3063  }
3064  }
3065  }
3066  }
3067  }
3068 
3069  else if( strcmp( save.chSave[ipPun],"LEIS") == 0)
3070  {
3071  if( strcmp(chTime,"LAST") != 0 )
3072  {
3073  /* get some column densities we shall need */
3074  double col_ci , col_oi , col_cii, col_heii;
3075  if( cdColm("carb" , 1 , &col_ci ) )
3076  TotalInsanity();
3077  if( cdColm("carb" , 2 , &col_cii ) )
3078  TotalInsanity();
3079  if( cdColm("oxyg" , 1 , &col_oi ) )
3080  TotalInsanity();
3081  if( cdColm("heli" , 2 , &col_heii ) )
3082  TotalInsanity();
3083  /* save Leiden structure - some numbers for the Leiden PDR model comparisons */
3084  fprintf( save.ipPnunit[ipPun],
3085  "%.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"
3086  "%.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"
3087  "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t"
3088  "%.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"
3089  "%.3e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
3090  /* depth of this point */
3091  radius.depth_mid_zone,
3092  /* A_V for an extended source */
3093  0.00,
3094  /* A_V for a point source */
3095  rfield.extin_mag_V_point,
3096  /* temperature */
3097  phycon.te ,
3098  dense.xIonDense[ipHYDROGEN][0],
3099  hmi.H2_total,
3100  dense.xIonDense[ipCARBON][0],
3101  dense.xIonDense[ipCARBON][1],
3102  dense.xIonDense[ipOXYGEN][0],
3103  findspecieslocal("CO")->den,
3104  findspecieslocal("O2")->den,
3105  findspecieslocal("CH")->den,
3106  findspecieslocal("OH")->den,
3107  dense.eden,
3108  dense.xIonDense[ipHELIUM][1],
3109  dense.xIonDense[ipHYDROGEN][1],
3110  findspecieslocal("H3+")->den,
3111  colden.colden[ipCOL_H0],
3112  colden.colden[ipCOL_H2g]+colden.colden[ipCOL_H2s],
3113  col_ci,
3114  col_cii,
3115  col_oi,
3116  findspecieslocal("CO")->column,
3117  findspecieslocal("O2")->column,
3118  findspecieslocal("CH")->column,
3119  findspecieslocal("OH")->column,
3120  colden.colden[ipCOL_elec],
3121  col_heii,
3122  colden.colden[ipCOL_Hp],
3123  colden.colden[ipCOL_H3p],
3124  hmi.H2_Solomon_dissoc_rate_used_H2g ,
3125  gv.rate_h2_form_grains_used_total,
3126  hmi.H2_photodissoc_used_H2g,
3127  hmi.UV_Cont_rel2_Draine_DB96_depth,
3128  /* CO and C dissociation rate */
3129  mole.findrk("PHOTON,CO=>C,O"),
3130  /* total CI ionization rate */
3131  ionbal.PhotoRate_Shell[ipCARBON][0][2][0],
3132  /* total heating, erg cm-3 s-1 */
3133  thermal.htot,
3134  /* total cooling, erg cm-3 s-1 */
3135  thermal.ctot,
3136  /* GrnP grain photo heating */
3137  thermal.heating[0][13],
3138  /* grain collisional cooling */
3139  MAX2(0.,gv.GasCoolColl),
3140  /* grain collisional heating */
3141  -1.*MIN2(0.,gv.GasCoolColl),
3142  /* COds - CO dissociation heating */
3143  thermal.heating[0][9],
3144  /* H2dH-Heating due to H2 dissociation */
3145  hmi.HeatH2Dish_used,
3146  /* H2vH-Heating due to collisions with H2 */
3147  hmi.HeatH2Dexc_used ,
3148  /* ChaT - charge transfer heating */
3149  thermal.heating[0][24] ,
3150  /* cosmic ray heating */
3151  thermal.heating[1][6] ,
3152  /* heating due to atoms of various heavy elements */
3153  thermal.heating[ipMAGNESIUM][0],
3154  thermal.heating[ipSULPHUR][0],
3155  thermal.heating[ipSILICON][0],
3156  thermal.heating[ipIRON][0],
3157  thermal.heating[ipSODIUM][0],
3158  thermal.heating[ipALUMINIUM][0],
3159  thermal.heating[ipCARBON][0],
3160  TauLines[ipT610].Coll().cool(),
3161  TauLines[ipT370].Coll().cool(),
3162  TauLines[ipT157].Coll().cool(),
3163  TauLines[ipT63].Coll().cool(),
3164  TauLines[ipT146].Coll().cool() );
3165  }
3166  }
3167 
3168 # ifdef USE_NLTE7
3169  else if( strcmp( save.chSave[ipPun],"NLTE") == 0)
3170  {
3171  if( strcmp(chTime,"LAST") == 0 )
3172  {
3173  //Generate output for NLTE7 conference
3174  runNLTE(ipPun);
3175  }
3176  }
3177 # endif
3178 
3179 
3180  /* Print out the FeII energy level file into Stout format*/
3181  else if( strcmp( save.chSave[ipPun],"LY1") == 0)
3182  {
3183  static bool runonce = true;
3184  if (runonce )
3185  {
3186  for( long ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
3187  {
3188  ASSERT(FeII.FeIINRGs[ipHi] >= 0.);
3189  ASSERT(FeII.FeIISTWT[ipHi] >= 0.);
3190  fprintf(save.ipPnunit[ipPun], "%li\t%10.3f\t%3.1f\n",ipHi+1,FeII.FeIINRGs[ipHi],FeII.FeIISTWT[ipHi]);
3191  }
3192  runonce = false;
3193  fprintf(save.ipPnunit[ipPun], "-1\n");
3194  }
3195  }
3196 
3197  /* Print out the FeII transition probablility file into Stout format*/
3198  else if( strcmp( save.chSave[ipPun],"LY2") == 0)
3199  {
3200  static bool runonce = true;
3201  if (runonce )
3202  {
3203  for( long ipLo = 0; ipLo < FeII.nFeIILevel_malloc; ipLo++)
3204  {
3205  for( long ipHi=0; ipHi < FeII.nFeIILevel_malloc; ipHi++ )
3206  {
3207  if( FeII.FeIIAul[ipHi][ipLo] >= 0. && FeII.FeIIAul[ipHi][ipLo] != 1e-5f && FeII.FeIIAul[ipHi][ipLo] != 1e-6f )
3208  {
3209  fprintf(save.ipPnunit[ipPun], "A\t%li\t%li\t%8.2e\n",
3210  ipLo+1,ipHi+1,FeII.FeIIAul[ipHi][ipLo]);
3211  }
3212  }
3213  }
3214  runonce = false;
3215  fprintf(save.ipPnunit[ipPun], "-1\n");
3216  }
3217  }
3218 
3219  /* Print out the FeII Collision data file into Stout format*/
3220  else if( strcmp( save.chSave[ipPun],"LY3") == 0)
3221  {
3222  static realnum tt[8]={1e3f,3e3f,5e3f,7e3f,1e4f,12e3f,15e3f,2e4f};
3223  static bool runonce = true;
3224  if (runonce )
3225  {
3226  fprintf(save.ipPnunit[ipPun], "TEMP\t8");
3227  for( int k = 0; k < 8; k++)
3228  {
3229  fprintf(save.ipPnunit[ipPun], "\t%8.2e",tt[k]);
3230  }
3231  fprintf(save.ipPnunit[ipPun], "\n");
3232  for( long ipHi = 1; ipHi < FeII.nFeIILevel_malloc; ipHi++)
3233  {
3234  for( long ipLo=0; ipLo < ipHi; ipLo++ )
3235  {
3236  bool skipLine = false;
3237  for( int k = 0; k < 8; k++)
3238  {
3239  if( FeII.FeIIColl[ipHi][ipLo][k] == -2. || FeII.FeIIColl[ipHi][ipLo][k] == 0.)
3240  {
3241  skipLine = true;
3242  break;
3243  }
3244  }
3245  if( !skipLine )
3246  {
3247  fprintf(save.ipPnunit[ipPun], "CSELECTRON\t%li\t%li",ipLo+1,ipHi+1);
3248  for( int k = 0; k < 8; k++)
3249  {
3250  fprintf(save.ipPnunit[ipPun], "\t%8.2e",FeII.FeIIColl[ipHi][ipLo][k]);
3251  }
3252  fprintf(save.ipPnunit[ipPun], "\n");
3253  }
3254  }
3255  }
3256  runonce = false;
3257  fprintf(save.ipPnunit[ipPun], "-1\n");
3258  }
3259  }
3260 
3261  else if( strcmp( save.chSave[ipPun],"LLST") == 0)
3262  {
3263  /* save linelist command - do on last iteration */
3264  if( strcmp(chTime,"LAST") == 0 )
3265  {
3266  fprintf( save.ipPnunit[ipPun], "iteration %li" , iteration );
3267  if( save.punarg[ipPun][1] )// column print
3268  fprintf( save.ipPnunit[ipPun], "\n" );
3269 
3270  /* -1 is flag saying that this save command was not set */
3271  if( save.nLineList[ipPun] < 0 )
3272  TotalInsanity();
3273 
3274  int LineType = 0;
3275  if( save.lgEmergent[ipPun] )
3276  LineType = 1;
3277  if( save.lgCumulative[ipPun] )
3278  LineType += 2;
3279 
3280  /* loop over all lines in the file we read */
3281  for( j=0; j<save.nLineList[ipPun]; ++j )
3282  {
3283  double relative , absolute, PrtQuantity;
3284  if( (cdLine( save.chLineListLabel[ipPun][j] ,
3285  save.wlLineList[ipPun][j] ,
3286  &relative , &absolute , LineType ) ) <=0 )
3287  {
3288  if( !h2.lgEnabled && strncmp( save.chLineListLabel[ipPun][j] , "H2 " , 4 )==0 )
3289  {
3290  static bool lgMustPrintFirstTime = true;
3291  if( lgMustPrintFirstTime )
3292  {
3293  /* it's an H2 line and H2 is not being done - ignore it */
3294  fprintf( ioQQQ,"Did not find an H2 line, the large model is not "
3295  "included, so I will ignore it. Log intensity set to -30.\n" );
3296  fprintf( ioQQQ,"I will totally ignore any future missed H2 lines\n");
3297  lgMustPrintFirstTime = false;
3298  }
3299  relative = -30.f;
3300  absolute = -30.f;
3301  }
3302  else if( lgAbort )
3303  {
3304  /* we are in abort mode */
3305  relative = -30.f;
3306  absolute = -30.f;
3307  }
3308  else
3309  {
3310  fprintf(ioQQQ,"DISASTER - did not find a line in the Line List table\n");
3311  cdEXIT(EXIT_FAILURE);
3312  }
3313  }
3314 
3315  /* options to do either relative or absolute intensity
3316  * default is relative, is absolute keyword on line then
3317  * punarg set to 1 */
3318  /* straight line intensities */
3319  if( save.punarg[ipPun][0] > 0 )
3320  PrtQuantity = pow(10. , absolute);
3321  else
3322  PrtQuantity = relative;
3323 
3324  // column mode, print label
3325  if( save.punarg[ipPun][1] )
3326  {
3327  /* if taking ratio then put div sign between pairs */
3328  if( save.lgLineListRatio[ipPun] && is_odd(j) )
3329  fprintf( save.ipPnunit[ipPun] , "/" );
3330 
3331  fprintf( save.ipPnunit[ipPun], "%s ", save.chLineListLabel[ipPun][j] );
3332  char chTemp[MAX_HEADER_SIZE];
3333  sprt_wl( chTemp, save.wlLineList[ipPun][j] );
3334  fprintf( save.ipPnunit[ipPun], "%s ", chTemp );
3335  }
3336 
3337  /* if taking ratio print every other line as ratio
3338  * with previous line */
3339  if( save.lgLineListRatio[ipPun] )
3340  {
3341  /* do line pair ratios */
3342  static double SaveQuantity = 0;
3343  if( is_odd(j) )
3344  fprintf( save.ipPnunit[ipPun], "\t%.4e" ,
3345  SaveQuantity / SDIV( PrtQuantity ) );
3346  else
3347  SaveQuantity = PrtQuantity;
3348  }
3349  else
3350  {
3351  fprintf( save.ipPnunit[ipPun], "\t%.4e" , PrtQuantity );
3352  }
3353  // column printout, but check if first of pair
3354  if( save.punarg[ipPun][1] )
3355  {
3356  if( !save.lgLineListRatio[ipPun] ||
3357  is_odd(j) )
3358  fprintf( save.ipPnunit[ipPun], "\n" );
3359  }
3360  }
3361  fprintf( save.ipPnunit[ipPun], "\n" );
3362  }
3363  }
3364 
3365  else if( strcmp( save.chSave[ipPun],"CHRT") == 0)
3366  {
3367  /* save chemistry rates command */
3368  if( strcmp(chTime,"LAST") != 0 )
3369  {
3370  bool lgHeader, lgData;
3371  if( save.lgPunHeader[ipPun] )
3372  {
3373  lgHeader = true;
3374  lgData = false;
3375  mole_punch(save.ipPnunit[ipPun],save.optname[ipPun].c_str(),save.chSaveArgs[ipPun],lgHeader,lgData,radius.depth_mid_zone);
3376  }
3377  save.lgPunHeader[ipPun] = false;
3378  lgHeader = false;
3379  lgData = true;
3380  mole_punch(save.ipPnunit[ipPun],save.optname[ipPun].c_str(),save.chSaveArgs[ipPun],lgHeader,lgData,radius.depth_mid_zone);
3381  }
3382  }
3383 
3384  else if( strcmp(save.chSave[ipPun],"MAP ") == 0 )
3385  {
3386  /* do the map now if we are at the zone, or if this
3387  * is the LAST call to this routine and map not done yet */
3388  if( !hcmap.lgMapDone &&
3389  (nzone == hcmap.MapZone || strcmp(chTime,"LAST") == 0 ) )
3390  {
3391  lgTlkSav = called.lgTalk;
3392  called.lgTalk = cpu.i().lgMPI_talk();
3393  hcmap.lgMapBeingDone = true;
3394  map_do(save.ipPnunit[ipPun] , " map");
3395  called.lgTalk = lgTlkSav;
3396  }
3397  }
3398 
3399  else if( strcmp(save.chSave[ipPun],"MOLE") == 0 )
3400  {
3401  if( save.lgPunHeader[ipPun] )
3402  {
3403  fprintf( save.ipPnunit[ipPun],
3404  "#molecular species will follow:\n");
3405  fprintf( save.ipPnunit[ipPun],
3406  "#depth\tAV(point)\tAV(extend)\tCO diss rate\tC recom rate");
3407 
3408  for(i=0; i<mole_global.num_calc; ++i )
3409  {
3410  fprintf( save.ipPnunit[ipPun], "\t%s", mole_global.list[i]->label.c_str() );
3411  }
3412  fprintf ( save.ipPnunit[ipPun], "\n");
3413  save.lgPunHeader[ipPun] = false;
3414  }
3415  if( strcmp(chTime,"LAST") != 0 )
3416  {
3417  /* molecules, especially for PDR, first give radius */
3418  fprintf( save.ipPnunit[ipPun], "%.5e\t" , radius.depth_mid_zone );
3419 
3420  /* visual extinction of point source (star)*/
3421  fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);
3422 
3423  /* visual extinction of an extended source (like a PDR)*/
3424  fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
3425 
3426  /* carbon monoxide photodissociation rate */
3427  fprintf( save.ipPnunit[ipPun], "%.5e\t" , mole.findrk("PHOTON,CO=>C,O") );
3428 
3429  /* carbon recombination rate */
3430  fprintf( save.ipPnunit[ipPun], "%.5e" , ionbal.RateRecomTot[ipCARBON][0] );
3431 
3432  /* now do all the molecules */
3433  for(j=0; j<mole_global.num_calc; ++j )
3434  {
3435  fprintf(save.ipPnunit[ipPun],"\t%.2e",mole.species[j].den );
3436  }
3437 
3438  fprintf(save.ipPnunit[ipPun],"\n");
3439  }
3440  }
3441 
3442  else if( strcmp(save.chSave[ipPun],"OPAC") == 0 )
3443  {
3444  /* save opacity- routine will parse which type of opacity save to do */
3445  if( save.lgSaveEveryZone[ipPun] || strcmp(chTime,"LAST") == 0 )
3446  save_opacity(save.ipPnunit[ipPun],ipPun);
3447  }
3448 
3449  /* save coarse optical depths command */
3450  else if( strcmp(save.chSave[ipPun],"OPTc") == 0 )
3451  {
3452  if( save.lgSaveEveryZone[ipPun] || strcmp(chTime,"LAST") == 0 )
3453  {
3454  for( j=0; j < rfield.nflux; j++ )
3455  {
3456  fprintf( save.ipPnunit[ipPun],
3457  "%13.5e\t%.3e\t%12.4e\t%.3e\n",
3458  AnuUnit(rfield.AnuOrg[j]),
3459  opac.TauAbsFace[j]+opac.TauScatFace[j],
3460  opac.TauAbsFace[j],
3461  opac.TauScatFace[j] );
3462  }
3463  }
3464  }
3465 
3466  /* save fine optical depths command */
3467  else if( strcmp(save.chSave[ipPun],"OPTf") == 0 )
3468  {
3469  if( save.lgSaveEveryZone[ipPun] || strcmp(chTime,"LAST") == 0 )
3470  {
3471  long nu_hi , nskip;
3472  if( save.punarg[ipPun][0] > 0. )
3473  /* possible lower bounds to energy range - will be zero if not set */
3474  j = ipFineCont( save.punarg[ipPun][0] );
3475  else
3476  j = 0;
3477 
3478  /* upper limit */
3479  if( save.punarg[ipPun][1]> 0. )
3480  nu_hi = ipFineCont( save.punarg[ipPun][1]);
3481  else
3482  nu_hi = rfield.nfine;
3483 
3484  /* we will bring nskip cells together into one printed
3485  * number to make output smaller - default is 10 */
3486  nskip = (long)save.punarg[ipPun][2];
3487  do
3488  {
3489  realnum sum1 = rfield.fine_opt_depth[j];
3490  realnum sum2 = rfield.fine_opac_zone[j];
3491  /* want to report the central wavelength of the cell */
3492  realnum xnu = rfield.fine_anu[j];
3493  for( jj=1; jj<nskip; ++jj )
3494  {
3495  sum1 += rfield.fine_opt_depth[j+jj];
3496  sum2 += rfield.fine_opac_zone[j+jj];
3497  xnu += rfield.fine_anu[j+jj];
3498  }
3499  if( sum2 > 0. )
3500  fprintf( save.ipPnunit[ipPun],
3501  "%12.6e\t%.3e\t%.3e\n",
3502  AnuUnit(xnu/nskip),
3503  sum1/nskip ,
3504  sum2/nskip);
3505  j += nskip;
3506  }while( j < nu_hi );
3507  }
3508  }
3509 
3510  else if( strcmp(save.chSave[ipPun]," OTS") == 0 )
3511  {
3512  ConMax = 0.;
3513  xLinMax = 0.;
3514  opConSum = 0.;
3515  opLinSum = 0.;
3516  ipLinMax = 1;
3517  ipConMax = 1;
3518 
3519  for( j=0; j < rfield.nflux; j++ )
3520  {
3521  opConSum += rfield.otscon[j]*opac.opacity_abs[j];
3522  opLinSum += rfield.otslin[j]*opac.opacity_abs[j];
3523  if( rfield.otslin[j]*opac.opacity_abs[j] > xLinMax )
3524  {
3525  xLinMax = rfield.otslin[j]*opac.opacity_abs[j];
3526  ipLinMax = j+1;
3527  }
3528  if( rfield.otscon[j]*opac.opacity_abs[j] > ConMax )
3529  {
3530  ConMax = rfield.otscon[j]*opac.opacity_abs[j];
3531  ipConMax = j+1;
3532  }
3533  }
3534  fprintf( save.ipPnunit[ipPun],
3535  "tot con lin=%.2e%.2e lin=%.4s%.4e%.2e con=%.4s%.4e%.2e\n",
3536  opConSum, opLinSum, rfield.chLineLabel[ipLinMax-1]
3537  , rfield.anu[ipLinMax-1], xLinMax, rfield.chContLabel[ipConMax-1]
3538  , rfield.anu[ipConMax-1], ConMax );
3539  }
3540 
3541  else if( strcmp(save.chSave[ipPun],"OVER") == 0 )
3542  {
3543  /* save overview
3544  * this is the floor for the smallest ionization fractions printed */
3545  double toosmall = SMALLFLOAT ,
3546  hold;
3547 
3548  /* overview of model results,
3549  * depth, te, hden, eden, ion fractions H, He, c, O */
3550  if( strcmp(chTime,"LAST") != 0 )
3551  {
3552 
3553  /* print the depth */
3554  fprintf( save.ipPnunit[ipPun], "%.5e\t", radius.depth_mid_zone );
3555 
3556  /* temperature, heating */
3557  if(dynamics.Cool() > dynamics.Heat())
3558  {
3559  fprintf( save.ipPnunit[ipPun], "%.4f\t%.3f",
3560  log10(phycon.te), log10(SDIV(thermal.htot-dynamics.Heat())) );
3561  }
3562  else
3563  {
3564  double diff = fabs(thermal.htot-dynamics.Cool());
3565  fprintf( save.ipPnunit[ipPun], "%.4f\t%.3f",
3566  log10(phycon.te), log10(SDIV(diff)) );
3567  }
3568 
3569  /* hydrogen and electron densities */
3570  fprintf( save.ipPnunit[ipPun], "\t%.4f\t%.4f",
3571  log10(dense.gas_phase[ipHYDROGEN]), log10(dense.eden) );
3572 
3573  /* molecular fraction of hydrogen */
3574  fprintf( save.ipPnunit[ipPun], "\t%.4f",
3575  /*log10(MAX2(toosmall,2.*findspecies("H2")->den/dense.gas_phase[ipHYDROGEN])) );*/
3576  log10(MAX2(toosmall,2.*hmi.H2_total/dense.gas_phase[ipHYDROGEN])) );
3577 
3578  /* ionization fractions of hydrogen */
3579  fprintf( save.ipPnunit[ipPun], "\t%.4f\t%.4f",
3580  log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN])),
3581  log10(MAX2(toosmall,dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN])) );
3582 
3583  /* ionization fractions of helium */
3584  for( j=1; j <= 3; j++ )
3585  {
3586  double arg1 = SDIV(dense.gas_phase[ipHELIUM]);
3587  arg1 = MAX2(toosmall,dense.xIonDense[ipHELIUM][j-1]/arg1 );
3588  fprintf( save.ipPnunit[ipPun], "\t%.4f",
3589  log10(arg1) );
3590  }
3591 
3592  /* carbon monoxide molecular fraction of CO */
3593  hold = SDIV(dense.gas_phase[ipCARBON]);
3594  hold = findspecieslocal("CO")->den/hold;
3595  hold = MAX2(toosmall, hold );
3596  fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(hold) );
3597 
3598  /* ionization fractions of carbon */
3599  for( j=1; j <= 4; j++ )
3600  {
3601  hold = SDIV(dense.gas_phase[ipCARBON]);
3602  hold = MAX2(toosmall,dense.xIonDense[ipCARBON][j-1]/hold);
3603  fprintf( save.ipPnunit[ipPun], "\t%.4f",
3604  log10(hold) );
3605  }
3606 
3607  /* ionization fractions of oxygen */
3608  for( j=1; j <= 6; j++ )
3609  {
3610  hold = SDIV(dense.gas_phase[ipOXYGEN]);
3611  hold = MAX2(toosmall,dense.xIonDense[ipOXYGEN][j-1]/hold);
3612  fprintf( save.ipPnunit[ipPun], "\t%.4f",
3613  log10(hold) );
3614  }
3615 
3616  // molecular fraction of H2O
3617  hold = SDIV(dense.gas_phase[ipOXYGEN]);
3618  hold = findspecieslocal("H2O")->den/hold;
3619  hold = MAX2(toosmall, hold );
3620  fprintf( save.ipPnunit[ipPun], "\t%.4f", log10(hold) );
3621 
3622  /* visual extinction of point source (star)*/
3623  fprintf( save.ipPnunit[ipPun], "\t%.2e" , rfield.extin_mag_V_point);
3624 
3625  /* visual extinction of an extended source (like a PDR)*/
3626  fprintf( save.ipPnunit[ipPun], "\t%.2e\n" , rfield.extin_mag_V_extended);
3627  }
3628  }
3629 
3630  else if( strcmp(save.chSave[ipPun]," PDR") == 0 )
3631  {
3632  /* this is the save PDR command */
3633  if( strcmp(chTime,"LAST") != 0 )
3634  {
3635  /* convert optical depth at wavelength of V filter
3636  * into magnitudes of extinction */
3637  /* >>chyng 03 feb 25, report extinction to illuminated face,
3638  * rather than total extinction which included far side when
3639  * sphere was set */
3640  /*av = opac.TauTotalGeo[0][rfield.ipV_filter-1]*1.08574;*/
3641 
3642  fprintf( save.ipPnunit[ipPun],
3643  "%.5e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t",
3644  radius.depth_mid_zone,
3645  /* total hydrogen column density, all forms */
3646  colden.colden[ipCOL_HTOT],
3647  phycon.te,
3648  /* fraction of H that is atomic */
3649  dense.xIonDense[ipHYDROGEN][0]/dense.gas_phase[ipHYDROGEN],
3650  /* ratio of n(H2) to total H, == 0.5 when fully molecular */
3651  2.*findspecieslocal("H2")->den/dense.gas_phase[ipHYDROGEN],
3652  2.*findspecieslocal("H2*")->den/dense.gas_phase[ipHYDROGEN],
3653  /* atomic to total carbon */
3654  dense.xIonDense[ipCARBON][0]/SDIV(dense.gas_phase[ipCARBON]),
3655  findspecieslocal("CO")->den/SDIV(dense.gas_phase[ipCARBON]),
3656  findspecieslocal("H2O")->den/SDIV(dense.gas_phase[ipOXYGEN]),
3657  /* hmi.UV_Cont_rel2_Habing_TH85 is field relative to Habing background, dimensionless */
3658  hmi.UV_Cont_rel2_Habing_TH85_depth);
3659 
3660  /* visual extinction due to dust alone, of point source (star)*/
3661  fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_point);
3662 
3663  /* visual extinction due to dust alone, of an extended source (like a PDR)*/
3664  fprintf( save.ipPnunit[ipPun], "%.2e\t" , rfield.extin_mag_V_extended);
3665 
3666  /* visual extinction (all sources) of a point source (like a PDR)*/
3667  fprintf( save.ipPnunit[ipPun], "%.2e\n" , opac.TauAbsGeo[0][rfield.ipV_filter] );
3668  }
3669  }
3670 
3671  /* performance characteristics per zone */
3672  else if( strcmp(save.chSave[ipPun],"PERF") == 0 )
3673  {
3674  if( strcmp(chTime,"LAST") != 0 )
3675  {
3676  static double ElapsedTime , ZoneTime;
3677  if( nzone<=1 )
3678  {
3679  ElapsedTime = cdExecTime();
3680  ZoneTime = 0.;
3681  }
3682  else
3683  {
3684  double t = cdExecTime();
3685  ZoneTime = t - ElapsedTime;
3686  ElapsedTime = t;
3687  }
3688 
3689  /* zone, time for this zone, elapsed time */
3690  fprintf( save.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\t%li",
3691  nzone, ZoneTime , ElapsedTime, conv.nPres2Ioniz );
3692  // print various loop counters
3693  for( long i=0; i<NTYPES; ++i )
3694  fprintf( save.ipPnunit[ipPun], "\t%li", conv.getCounterZone(i) );
3695  fprintf( save.ipPnunit[ipPun], "\n" );
3696  }
3697  }
3698 
3699  else if( strcmp(save.chSave[ipPun],"PHYS") == 0 )
3700  {
3701  if( strcmp(chTime,"LAST") != 0 )
3702  {
3703  /* save physical conditions */
3704  fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.3e\t%.3e\t%.3e\t%.3e\t%.3e\n",
3705  radius.depth_mid_zone, phycon.te, dense.gas_phase[ipHYDROGEN],
3706  dense.eden, thermal.htot, wind.AccelTotalOutward, geometry.FillFac );
3707  }
3708  }
3709 
3710  else if( strcmp(save.chSave[ipPun],"PRES") == 0 )
3711  {
3712  /* the save pressure command */
3713  if( strcmp(chTime,"LAST") != 0 )
3714  {
3715  fprintf( save.ipPnunit[ipPun],
3716  "%.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",
3717  /*A 1 #P depth */
3718  radius.depth_mid_zone,
3719  /*B 2 Perror */
3720  pressure.PresTotlError*100.,
3721  /*C 3 Pcurrent */
3722  pressure.PresTotlCurr,
3723  /*D 4 Pln + pintg
3724  * >>chng 06 apr 19, subtract pinzon the acceleration added in this zone
3725  * since is not total at outer edge of zone, above is at inner edge */
3726  pressure.PresTotlInit + pressure.PresInteg - pressure.pinzon,
3727  /*E 5 pgas (0) */
3728  pressure.PresTotlInit,
3729  /*F 6 Pgas */
3730  pressure.PresGasCurr,
3731  /*G 7 Pram */
3732  pressure.PresRamCurr,
3733  /*H 8 P rad in lines */
3734  pressure.pres_radiation_lines_curr,
3735  /*I 9 Pinteg subtract continuum rad pres which has already been added on */
3736  pressure.PresInteg - pressure.pinzon,
3737  /*J 10 V(wind km/s) wind speed in km/s */
3738  wind.windv/1e5,
3739  /*K cad(km/s) sound speed in km/s */
3740  timesc.sound_speed_adiabatic/1e5,
3741  /* the magnetic pressure */
3742  magnetic.pressure ,
3743  /* the local turbulent velocity in km/s */
3744  DoppVel.TurbVel/1e5 ,
3745  /* turbulent pressure */
3746  pressure.PresTurbCurr*DoppVel.lgTurb_pressure ,
3747  /* gravitational pressure */
3748  pressure.IntegRhoGravity,
3749  // the integral of electron scattering acceleration in
3750  // the absence of any absorptio, minus acceleration in current
3751  // zone, which has been added in - done this way, result is
3752  // zero in first zonen
3753  pressure.PresIntegElecThin-pressure.pinzon_PresIntegElecThin,
3754  // is this converged?
3755  TorF(conv.lgConvPres) );
3756  }
3757  }
3758  else if( strcmp(save.chSave[ipPun],"PREL") == 0 )
3759  {
3760  /* line pressure contributors */
3761  fprintf( save.ipPnunit[ipPun],
3762  "%.5e\t%.3e\t%.3e\t",
3763  /*A 1 #P depth */
3764  radius.depth_mid_zone ,
3765  pressure.PresTotlCurr,
3766  pressure.pres_radiation_lines_curr/SDIV(pressure.PresTotlCurr) );
3767  PrtLinePres(save.ipPnunit[ipPun]);
3768 
3769  }
3770 
3771  else if( save.chSave[ipPun][0]=='R' )
3772  {
3773  /* work around internal limits to Microsoft vs compiler */
3774  if( strcmp(save.chSave[ipPun],"RADI") == 0 )
3775  {
3776  /* save radius information for all zones */
3777  if( strcmp(chTime,"LAST") != 0 )
3778  {
3779  fprintf( save.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n",
3780  nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
3781  radius.drad );
3782  }
3783  }
3784 
3785  else if( strcmp(save.chSave[ipPun],"RADO") == 0 )
3786  {
3787  /* save radius information for only the last zone */
3788  if( strcmp(chTime,"LAST") == 0 )
3789  {
3790  fprintf( save.ipPnunit[ipPun], "%ld\t%.5e\t%.4e\t%.4e\n",
3791  nzone, radius.Radius_mid_zone, radius.depth_mid_zone,
3792  radius.drad );
3793  }
3794  }
3795 
3796  else if( strcmp(save.chSave[ipPun],"RESU") == 0 )
3797  {
3798  /* save results of the calculation */
3799  if( strcmp(chTime,"LAST") == 0 )
3800  SaveResults(save.ipPnunit[ipPun]);
3801  }
3802  else
3803  {
3804  /* this can't happen */
3805  TotalInsanity();
3806  }
3807  }
3808 
3809  else if( strcmp(save.chSave[ipPun],"SECO") == 0 )
3810  {
3811  /* save secondary ionization */
3812  if( strcmp(chTime,"LAST") != 0 )
3813  fprintf(save.ipPnunit[ipPun],
3814  "%.5e\t%.3e\t%.3e\t%.3e\n",
3815  radius.depth ,
3816  secondaries.csupra[ipHYDROGEN][0],
3817  secondaries.csupra[ipHYDROGEN][0]*2.02,
3818  secondaries.x12tot );
3819  }
3820 
3821  else if( strcmp(save.chSave[ipPun],"SOUS") == 0 )
3822  {
3823  /* full spectrum of continuum source function at 1 depth
3824  * command was "save source spectrum" */
3825  if( strcmp(chTime,"LAST") != 0 )
3826  {
3827  limit = MIN2(rfield.ipMaxBolt,rfield.nflux);
3828  for( j=0; j < limit; j++ )
3829  {
3830  fprintf( save.ipPnunit[ipPun],
3831  "%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
3832  AnuUnit(rfield.AnuOrg[j]),
3833  rfield.ConEmitLocal[nzone][j]/rfield.widflx[j],
3834  opac.opacity_abs[j],
3835  rfield.ConSourceFcnLocal[nzone][j],
3836  rfield.ConSourceFcnLocal[nzone][j]/plankf(j) ,
3837  safe_div(rfield.ConSourceFcnLocal[nzone][j],rfield.flux[0][j]));
3838  }
3839  }
3840  }
3841 
3842  else if( strcmp(save.chSave[ipPun],"SOUD") == 0 )
3843  {
3844  /* parts of continuum source function vs depth
3845  * command was save source function depth */
3846  j = iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon + 2;
3847  fprintf( save.ipPnunit[ipPun],
3848  "%.4e\t%.4e\t%.4e\t%.4e\n",
3849  opac.TauAbsFace[j-1],
3850  rfield.ConEmitLocal[nzone][j-1]/rfield.widflx[j-1]/MAX2(1e-35,opac.opacity_abs[j-1]),
3851  rfield.otscon[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1],
3852  rfield.otscon[iso_sp[ipH_LIKE][ipHYDROGEN].fb[0].ipIsoLevNIonCon-1]/opac.opacity_abs[iso_sp[ipH_LIKE][ipHYDROGEN].fb[ipH1s].ipIsoLevNIonCon-1] );
3853  }
3854 
3855  /* this is save special option */
3856  else if( strcmp(save.chSave[ipPun],"SPEC") == 0 )
3857  {
3858  SaveSpecial(save.ipPnunit[ipPun],chTime);
3859  }
3860 
3861  /* this is save species option */
3862  else if( strcmp(save.chSave[ipPun],"SPCS") == 0 )
3863  {
3864  if( ( strcmp(chTime,"LAST") != 0 && strcmp(save.chSaveArgs[ipPun],"COLU") != 0 ) ||
3865  ( strcmp(chTime,"LAST") == 0 && strcmp(save.chSaveArgs[ipPun],"COLU") == 0 ) )
3866  SaveSpecies(save.ipPnunit[ipPun] , ipPun);
3867  }
3868 
3869  else if( strcmp(save.chSave[ipPun],"TEMP") == 0 )
3870  {
3871  static double deriv_old=-1;
3872  double deriv=-1. , deriv_sec;
3873  /* temperature and its derivatives */
3874  fprintf( save.ipPnunit[ipPun], "%.5e\t%.4e\t%.2e",
3875  radius.depth_mid_zone,
3876  phycon.te,
3877  thermal.dCooldT );
3878  /* if second zone then have one deriv */
3879  if( nzone >1 )
3880  {
3881  deriv = (phycon.te - struc.testr[nzone-2])/ radius.drad;
3882  fprintf( save.ipPnunit[ipPun], "\t%.2e", deriv );
3883  /* if third zone then have second deriv */
3884  if( nzone > 2 )
3885  {
3886  deriv_sec = (deriv-deriv_old)/ radius.drad;
3887  fprintf( save.ipPnunit[ipPun], "\t%.2e",
3888  deriv_sec );
3889  }
3890  deriv_old = deriv;
3891  }
3892  fprintf( save.ipPnunit[ipPun], "\n");
3893  }
3894 
3895  /* time dependent model */
3896  else if( strcmp(save.chSave[ipPun],"TIMD") == 0 )
3897  {
3898  if( strcmp(chTime,"LAST") == 0 )
3899  DynaPunchTimeDep( save.ipPnunit[ipPun] , "END" );
3900  }
3901 
3902  /* execution time per zone */
3903  else if( strcmp(save.chSave[ipPun],"XTIM") == 0 )
3904  {
3905  static double ElapsedTime , ZoneTime;
3906  if( nzone<=1 )
3907  {
3908  ElapsedTime = cdExecTime();
3909  ZoneTime = 0.;
3910  }
3911  else
3912  {
3913  double t = cdExecTime();
3914  ZoneTime = t - ElapsedTime;
3915  ElapsedTime = t;
3916  }
3917 
3918  /* zone, time for this zone, elapsed time */
3919  fprintf( save.ipPnunit[ipPun], " %ld\t%.3f\t%.2f\n",
3920  nzone, ZoneTime , ElapsedTime );
3921  }
3922 
3923  else if( strcmp(save.chSave[ipPun],"TPRE") == 0 )
3924  {
3925  /* temperature and its predictors, turned on with save tprid */
3926  fprintf( save.ipPnunit[ipPun], "%5ld %11.4e %11.4e %11.4e %g\n",
3927  nzone, phycon.TeInit, phycon.TeProp, phycon.te,
3928  (phycon.TeProp- phycon.te)/phycon.te );
3929  }
3930 
3931  else if( strcmp(save.chSave[ipPun],"WIND") == 0 )
3932  {
3933  /* wind velocity, radiative acceleration, and ratio total
3934  * to electron scattering acceleration */
3935  /* first test only save last zone */
3936  if( (save.punarg[ipPun][0] == 0 && strcmp(chTime,"LAST") == 0)
3937  ||
3938  /* this test save all zones */
3939  (save.punarg[ipPun][0] == 1 && strcmp(chTime,"LAST") != 0 ) )
3940  {
3941  fprintf( save.ipPnunit[ipPun],
3942  "%.5e\t%.5e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
3943  radius.Radius_mid_zone,
3944  radius.depth_mid_zone,
3945  wind.windv,
3946  wind.AccelTotalOutward,
3947  wind.AccelLine,
3948  wind.AccelCont ,
3949  wind.fmul ,
3950  wind.AccelGravity );
3951  }
3952  }
3953 
3954  else if( strcmp(save.chSave[ipPun],"XATT") == 0 )
3955  {
3956  /* attenuated incident continuum */
3957  ASSERT( grid.lgOutputTypeOn[2] );
3958 
3959  if( strcmp(chTime,"LAST") == 0 )
3960  {
3961  if( grid.lgGrid )
3962  saveFITSfile( save.ipPnunit[ipPun], 2 );
3963  else
3964  {
3965  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
3966  cdEXIT(EXIT_FAILURE);
3967  }
3968  }
3969  }
3970  else if( strcmp(save.chSave[ipPun],"XRFI") == 0 )
3971  {
3972  /* reflected incident continuum */
3973  ASSERT( grid.lgOutputTypeOn[3] );
3974 
3975  if( strcmp(chTime,"LAST") == 0 )
3976  {
3977  if( grid.lgGrid )
3978  saveFITSfile( save.ipPnunit[ipPun], 3 );
3979  else
3980  {
3981  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
3982  cdEXIT(EXIT_FAILURE);
3983  }
3984  }
3985  }
3986  else if( strcmp(save.chSave[ipPun],"XINC") == 0 )
3987  {
3988  /* incident continuum */
3989  ASSERT( grid.lgOutputTypeOn[1] );
3990 
3991  if( strcmp(chTime,"LAST") == 0 )
3992  {
3993  if( grid.lgGrid )
3994  saveFITSfile( save.ipPnunit[ipPun], 1 );
3995  else
3996  {
3997  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
3998  cdEXIT(EXIT_FAILURE);
3999  }
4000  }
4001  }
4002  else if( strcmp(save.chSave[ipPun],"XDFR") == 0 )
4003  {
4004  /* reflected diffuse continuous emission */
4005  ASSERT( grid.lgOutputTypeOn[5] );
4006 
4007  if( strcmp(chTime,"LAST") == 0 )
4008  {
4009  if( grid.lgGrid )
4010  saveFITSfile( save.ipPnunit[ipPun], 5 );
4011  else
4012  {
4013  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4014  cdEXIT(EXIT_FAILURE);
4015  }
4016  }
4017  }
4018  else if( strcmp(save.chSave[ipPun],"XDFO") == 0 )
4019  {
4020  /* diffuse continuous emission outward */
4021  ASSERT( grid.lgOutputTypeOn[4] );
4022 
4023  if( strcmp(chTime,"LAST") == 0 )
4024  {
4025  if( grid.lgGrid )
4026  saveFITSfile( save.ipPnunit[ipPun], 4 );
4027  else
4028  {
4029  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4030  cdEXIT(EXIT_FAILURE);
4031  }
4032  }
4033  }
4034  else if( strcmp(save.chSave[ipPun],"XLNR") == 0 )
4035  {
4036  /* reflected lines */
4037  ASSERT( grid.lgOutputTypeOn[7] );
4038 
4039  if( strcmp(chTime,"LAST") == 0 )
4040  {
4041  if( grid.lgGrid )
4042  saveFITSfile( save.ipPnunit[ipPun], 7 );
4043  else
4044  {
4045  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4046  cdEXIT(EXIT_FAILURE);
4047  }
4048  }
4049  }
4050  else if( strcmp(save.chSave[ipPun],"XLNO") == 0 )
4051  {
4052  /* outward lines */
4053  ASSERT( grid.lgOutputTypeOn[6] );
4054 
4055  if( strcmp(chTime,"LAST") == 0 )
4056  {
4057  if( grid.lgGrid )
4058  saveFITSfile( save.ipPnunit[ipPun], 6 );
4059  else
4060  {
4061  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4062  cdEXIT(EXIT_FAILURE);
4063  }
4064  }
4065  }
4066  else if( strcmp(save.chSave[ipPun],"XREF") == 0 )
4067  {
4068  /* total reflected, lines and continuum */
4069  ASSERT( grid.lgOutputTypeOn[9] );
4070 
4071  if( strcmp(chTime,"LAST") == 0 )
4072  {
4073  if( grid.lgGrid )
4074  saveFITSfile( save.ipPnunit[ipPun], 9 );
4075  else
4076  {
4077  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4078  cdEXIT(EXIT_FAILURE);
4079  }
4080  }
4081  }
4082  else if( strcmp(save.chSave[ipPun],"XTOT") == 0 )
4083  {
4084  /* total spectrum, reflected plus transmitted */
4085  ASSERT( grid.lgOutputTypeOn[0] );
4086 
4087  if( strcmp(chTime,"LAST") == 0 )
4088  {
4089  if( grid.lgGrid )
4090  saveFITSfile( save.ipPnunit[ipPun], 0 );
4091  else
4092  {
4093  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4094  cdEXIT(EXIT_FAILURE);
4095  }
4096  }
4097  }
4098  else if( strcmp(save.chSave[ipPun],"XTRN") == 0 )
4099  {
4100  /* total outward, lines and continuum */
4101  ASSERT( grid.lgOutputTypeOn[8] );
4102 
4103  if( strcmp(chTime,"LAST") == 0 )
4104  {
4105  if( grid.lgGrid )
4106  saveFITSfile( save.ipPnunit[ipPun], 8 );
4107  else
4108  {
4109  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4110  cdEXIT(EXIT_FAILURE);
4111  }
4112  }
4113  }
4114  else if( strcmp(save.chSave[ipPun],"XSPM") == 0 )
4115  {
4116  /* exp(-tau) to the illuminated face */
4117  ASSERT( grid.lgOutputTypeOn[10] );
4118 
4119  if( strcmp(chTime,"LAST") == 0 )
4120  {
4121  if( grid.lgGrid )
4122  saveFITSfile( save.ipPnunit[ipPun], 10 );
4123  else
4124  {
4125  fprintf( ioQQQ," Cannot save xspec files unless doing a grid.\n" );
4126  cdEXIT(EXIT_FAILURE);
4127  }
4128  }
4129  }
4130  // termination of second set of nested if's
4131  // error if we have not matched key
4132  /* there are a few "save" commands that are handled elsewhere
4133  * save dr is an example. These will have lgRealSave set false */
4134  // lgNoHitFirstBranch says did not find in previous nest of if's
4135  else if( save.lgRealSave[ipPun] && lgNoHitFirstBranch )
4136  {
4137  /* this is insanity, internal flag set in ParseSave not seen here */
4138  fprintf( ioQQQ, " PROBLEM DISASTER SaveDo does not recognize flag %4.4s set by ParseSave. This is impossible.\n",
4139  save.chSave[ipPun] );
4140  TotalInsanity();
4141  }
4142 
4143  /* print special hash string to separate out various iterations
4144  * chTime is LAST on last iteration
4145  * save.lgHashEndIter flag is true by default, set false
4146  * with "no hash" keyword on save command
4147  * save.lg_separate_iterations is true by default, set false
4148  * when save time dependent calcs since do not want special
4149  * character between time steps
4150  * grid.lgGrid is only true when doing a grid of calculations */
4151  if( strcmp(chTime,"LAST") == 0 &&
4152  !(iterations.lgLastIt && !grid.lgGrid ) &&
4153  save.lgHashEndIter[ipPun] &&
4154  save.lg_separate_iterations[ipPun] &&
4155  !save.lgFITS[ipPun] )
4156  {
4157  if( dynamics.lgTimeDependentStatic && strcmp( save.chHashString , "TIME_DEP" )==0 )
4158  {
4159  fprintf( save.ipPnunit[ipPun], "\"time=%f\n",
4160  dynamics.time_elapsed );
4161  }
4162  else
4163  {
4164  fprintf( save.ipPnunit[ipPun], "%s",
4165  save.chHashString );
4166  if( grid.lgGrid && ( iterations.lgLastIt || lgAbort ) )
4167  fprintf( save.ipPnunit[ipPun], " GRID_DELIMIT -- grid%09ld",
4168  optimize.nOptimiz );
4169  fprintf( save.ipPnunit[ipPun], "\n" );
4170  }
4171  }
4172  if( save.lgFLUSH )
4173  fflush( save.ipPnunit[ipPun] );
4174  }
4175  }
4176  return;
4177 }
4178 
4179 /*SaveLineIntensity produce the 'save lines intensity' output */
4180 STATIC void SaveLineIntensity(FILE * ioPUN, long int ipPun , realnum Threshold )
4181 {
4182  long int i;
4183 
4184  DEBUG_ENTRY( "SaveLineIntensity()" );
4185 
4186  /* used to save out all the emission line intensities
4187  * first initialize the line image reader */
4188 
4189  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4190  input.echo(ioPUN);
4191 
4192  /* now print any cautions or warnings */
4193  cdWarnings( ioPUN);
4194  cdCautions( ioPUN);
4195  fprintf( ioPUN, "zone=%5ld\n", nzone );
4196  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4197  fprintf( ioPUN, "begin emission lines\n" );
4198 
4199  /* only save non-zero intensities */
4200  SaveResults1Line(ioPUN," ",0,0.,"Start");
4201 
4202  // check whether intrinsic or emergent line emissivity
4203  bool lgEmergent = false;
4204  if( save.punarg[ipPun][0] > 0 )
4205  lgEmergent = true;
4206 
4207  for( i=0; i < LineSave.nsum; i++ )
4208  {
4209  // Threshold is zero by default on save line intensity,
4210  // all option sets to negative number so that we report all lines
4211  if( LineSv[i].SumLine[lgEmergent] > Threshold )
4212  {
4213  SaveResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength,
4214  LineSv[i].SumLine[save.lgEmergent[ipPun]], "Line ");
4215  }
4216  }
4217 
4218  SaveResults1Line(ioPUN," ",0,0.,"Flush");
4219 
4220  fprintf( ioPUN, " \n" );
4221  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4222 
4223  return;
4224 }
4225 
4226 /* lgSaveOpticalDepths true says save optical depths */
4227 static bool lgPopsFirstCall , lgSaveOpticalDepths;
4228 
4229 /*SaveLineStuff save optical depths or source functions for all transferred lines */
4230 STATIC void SaveLineStuff(
4231  FILE * ioPUN,
4232  const char *chJob ,
4233  realnum xLimit )
4234 {
4235 
4236  long int nelem , ipLo , ipHi , i , ipISO;
4237  long index = 0;
4238  static bool lgFirst=true;
4239 
4240  DEBUG_ENTRY( "SaveLineStuff()" );
4241 
4242  /*find out which job this is and set a flag to use later */
4243  if( strcmp( &*chJob , "optical" ) == 0 )
4244  {
4245  /* save line optical depths */
4246  lgSaveOpticalDepths = true;
4247  lgPopsFirstCall = false;
4248  }
4249  else if( strcmp( &*chJob , "populat" ) == 0 )
4250  {
4251  lgSaveOpticalDepths = false;
4252  /* level population information */
4253  if( lgFirst )
4254  {
4255  lgPopsFirstCall = true;
4256  fprintf(ioPUN,"index\tAn.ion\tgLo\tgUp\tE(wn)\tgf\n");
4257  lgFirst = false;
4258  }
4259  else
4260  {
4261  lgPopsFirstCall = false;
4262  }
4263  }
4264  else
4265  {
4266  fprintf( ioQQQ, " insane job in SaveLineStuff =%s\n",
4267  &*chJob );
4268  cdEXIT(EXIT_FAILURE);
4269  }
4270 
4271  /* loop over all lines, calling put1Line to create info (routine located below) */
4272  /* hydrogen like lines */
4273  /* >>chng 02 may 16, had been explicit H and He-like loops */
4274  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
4275  {
4276  for( nelem=ipISO; nelem < LIMELM; nelem++ )
4277  {
4278  if( dense.lgElmtOn[nelem] )
4279  {
4280  /* 06 aug 28, from numLevels_max to _local. */
4281  for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
4282  {
4283  for( ipLo=0; ipLo <ipHi; ipLo++ )
4284  {
4285  if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
4286  continue;
4287 
4288  ++index;
4289  Save1Line( iso_sp[ipISO][nelem].trans(ipHi,ipLo), ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) );
4290  }
4291  }
4292  /* also do extra Lyman lines if optical depths are to be done,
4293  * these are line that are included only for absorption, not in the
4294  * model atoms */
4295  if( lgSaveOpticalDepths )
4296  {
4297  /* >>chng 02 aug 23, for he-like, had starting on much too high a level since
4298  * index was number of levels - caught by Adrian Turner */
4299  /* now output extra line lines, starting one above those already done above */
4300  /*for( ipHi=iso_sp[ipISO][nelem].numLevels_max; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )*/
4301  /* 06 aug 28, from numLevels_max to _local. */
4302  for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
4303  {
4304  ++index;
4305  Save1Line( ExtraLymanLines[ipISO][nelem][ipExtraLymanLines[ipISO][nelem][ipHi]], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[nelem]) );
4306  }
4307  }
4308  }
4309  }
4310  }
4311 
4312  /* index from 1 to skip over dummy line */
4313  for( i=1; i < nLevel1; i++ )
4314  {
4315  ++index;
4316  Save1Line( TauLines[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
4317  }
4318 
4319  for( i=0; i < nWindLine; i++ )
4320  {
4321  if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
4322  {
4323  ++index;
4324  Save1Line( TauLine2[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
4325  }
4326  }
4327 
4328  for( i=0; i < nUTA; i++ )
4329  {
4330  ++index;
4331  Save1Line( UTALines[i], ioPUN, xLimit, index, GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) );
4332  }
4333 
4334 
4335  /* do optical depths of FeII lines, if large atom is turned on */
4336  FeIIPunchLineStuff( ioPUN , xLimit , index);
4337 
4338  /* save optical depths of H2 lines */
4339  h2.H2_PunchLineStuff( ioPUN , xLimit , index);
4340 
4341  /*fprintf(ioPUN, "##################################\n"); */
4342  fprintf( ioPUN , "%s\n",save.chHashString );
4343  return;
4344 }
4345 
4346 /*Save1Line called by SaveLineStuff to produce output for one line */
4347 void Save1Line( const TransitionProxy& t , FILE * ioPUN , realnum xLimit , long index, realnum DopplerWidth )
4348 {
4349 
4350  if( lgSaveOpticalDepths )
4351  {
4352  /* optical depths, no special first time, only print them */
4353  if( t.Emis().TauIn() >= xLimit )
4354  {
4355  /* label like "C 4" or "H 1" */
4356  fprintf( ioPUN , "%2.2s%2.2s\t",
4357  elementnames.chElementSym[(*t.Hi()).nelem()-1] ,
4358  elementnames.chIonStage[(*t.Hi()).IonStg()-1] );
4359 
4360  /* print wavelengths, either line in main printout labels,
4361  * or in various units in exponential notation - prt_wl is in prt.c */
4362  if( strcmp( save.chConPunEnr[save.ipConPun], "labl" )== 0 )
4363  {
4364  prt_wl( ioPUN , t.WLAng() );
4365  }
4366  else
4367  {
4368  /* this converts energy in Rydbergs into any of the other units */
4369  fprintf( ioPUN , "%.7e", AnuUnit((realnum)(t.EnergyRyd())) );
4370  }
4371  /* print the optical depth */
4372  fprintf( ioPUN , "\t%.3f", t.Emis().TauIn() );
4373  /* damping constant */
4374  fprintf(ioPUN, "\t%.3e",
4375  t.Emis().dampXvel() / DopplerWidth );
4376  fprintf(ioPUN, "\n");
4377  }
4378  }
4379  else if( lgPopsFirstCall )
4380  {
4381  char chLbl[11];
4382  /* first call to line populations, print atomic parameters and indices */
4383  strcpy( chLbl, chLineLbl(t) );
4384  fprintf(ioPUN, "%li\t%s" , index , chLbl );
4385  /* stat weights */
4386  fprintf(ioPUN, "\t%.0f\t%.0f",
4387  (*t.Lo()).g() ,(*t.Hi()).g());
4388  /* energy difference, gf */
4389  fprintf(ioPUN, "\t%.2f\t%.3e",
4390  t.EnergyWN() ,t.Emis().gf());
4391  fprintf(ioPUN, "\n");
4392  }
4393  else
4394  {
4395  /* not first call, so do level populations and indices defined above */
4396  if( (*t.Hi()).Pop() > xLimit )
4397  {
4398  /* >>chng 05 may 08, add abundances, which for iso-seq species is
4399  * the density of the parent ion, for other lines, is unity.
4400  * had not been included so pops for iso seq were rel to parent ion.
4401  * caught by John Everett */
4402  /* multiplication by abundance no longer necessary since iso pops now denormalized */
4403  fprintf(ioPUN,"%li\t%.2e\t%.2e\n", index, (*t.Lo()).Pop(), (*t.Hi()).Pop() );
4404  }
4405  }
4406 }
4407 
4408 /*SaveNewContinuum produce the 'save new continuum' output */
4409 STATIC void SaveNewContinuum(FILE * ioPUN )
4410 {
4411  long int ipLo, ipHi,
4412  j ,
4413  ncells;
4414 
4415  double wllo, wlhi;
4416 
4417  double *energy,
4418  *cont_incid,
4419  *cont_atten,
4420  *diffuse_in,
4421  *diffuse_out,
4422  *emis_lines_out,
4423  *emis_lines_in;
4424 
4425  /* get the low limit */
4426  wllo = rfield.anu[0];
4427  /* get high-energy limit */
4428  wlhi = rfield.anu[rfield.nflux-1];
4429  /* use native continuum mesh */
4430  ipLo = ipoint(wllo)-1;
4431  ipHi = ipoint(wlhi)-1;
4432  ncells = ipHi - ipLo + 1;
4433 
4434  /* now allocate the space */
4435  energy = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4436  cont_incid = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4437  cont_atten = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4438  diffuse_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4439  diffuse_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4440  emis_lines_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4441  emis_lines_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double) );
4442  /*emis_lines_pump_out = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));
4443  emis_lines_pump_in = (double *)MALLOC( (size_t)(ncells+1)*sizeof(double));*/
4444 
4445  /* now convert to rydbergs */
4446  for( j=0; j<ncells; ++j )
4447  {
4448  energy[j] = rfield.AnuOrg[j+ipLo] - rfield.widflx[j+ipLo]/2.;
4449  }
4450 
4451  fixit();
4452  /* all of these should be rerouted to cdSPEC2, but units not right
4453  * need ergs multiplied for one, and continuum and lines may not be added correctly.
4454  * Goal is to abandon current cdSPEC and replace it with current cdSPEC2 */
4455 
4456 #if 1
4457  /* for cdSPEC the energy vector is the lower edge of the energy cell */
4458  /* get incident continuum */
4459  cdSPEC( 1 , ncells , cont_incid );
4460  /* get attenuated incident continuum */
4461  cdSPEC( 2 , ncells , cont_atten );
4462  /* get diffuse continuous emission, reflected */
4463  cdSPEC( 5 , ncells , diffuse_in );
4464  /* get continuous emission outward direction */
4465  cdSPEC( 4 , ncells , diffuse_out );
4467  /* get all outward lines */
4468  cdSPEC( 6 , ncells , emis_lines_out );
4469  /* get all reflected lines */
4470  cdSPEC( 7 , ncells , emis_lines_in );
4471 #else
4472  cdSPEC2( 1, rfield.nflux, 0, rfield.nflux - 1, cont_incid );
4473  /* get attenuated incident continuum */
4474  cdSPEC2( 2, rfield.nflux, 0, rfield.nflux - 1, cont_atten );
4475  /* get diffuse continuous emission, reflected */
4476  cdSPEC2( 5, rfield.nflux, 0, rfield.nflux - 1, diffuse_in );
4477  /* get continuous emission outward direction */
4478  cdSPEC2( 4, rfield.nflux, 0, rfield.nflux - 1, diffuse_out );
4479  /* get all outward lines */
4480  cdSPEC2( 6, rfield.nflux, 0, rfield.nflux - 1, emis_lines_out );
4481  /* get all reflected lines */
4482  cdSPEC2( 7, rfield.nflux, 0, rfield.nflux - 1, emis_lines_in );
4483 #endif
4484 
4485  /* for this example we will do a wavelength range */
4486  for( j=0; j<ncells-1; ++j )
4487  {
4488  /* photon energy in appropriate energy or wavelength units */
4489  fprintf( ioPUN,"%.5e\t", AnuUnit((realnum)(energy[j]+rfield.widflx[j+ipLo]/2.) ) );
4490  fprintf( ioPUN,"%.3e\t", cont_incid[j] );
4491  fprintf( ioPUN,"%.3e\t", cont_atten[j] );
4492  fprintf( ioPUN,"%.3e\t", diffuse_in[j]+diffuse_out[j] );
4493  fprintf( ioPUN,"%.3e",
4494  emis_lines_out[j]+emis_lines_in[j]/*+emis_lines_pump_out[j]+emis_lines_pump_in[j]*/ );
4495  fprintf( ioPUN,"\n" );
4496  }
4497 
4498  free(energy);
4499  free(cont_incid);
4500  free(diffuse_in);
4501  free(diffuse_out);
4502  free(cont_atten);
4503  free(emis_lines_out);
4504  free(emis_lines_in);
4505  /*free(emis_lines_pump_out);
4506  free(emis_lines_pump_in );*/
4507 }
4508 
4509 /* save AGN3 hemiss, for Chapter 4, routine is below */
4510 STATIC void AGN_Hemis(FILE *ioPUN )
4511 {
4512  const int NTE = 4;
4513  realnum te[NTE] = {5000., 10000., 15000., 20000.};
4514  realnum *agn_continuum[NTE];
4515  double TempSave = phycon.te;
4516  long i , j;
4517 
4518  DEBUG_ENTRY( "AGN_Hemis()" );
4519 
4520  /* make table of continuous emission at various temperatuers */
4521  /* first allocate space */
4522  for( i=0;i<NTE; ++i)
4523  {
4524  agn_continuum[i] = (realnum *)MALLOC((unsigned)rfield.nflux*sizeof(realnum) );
4525 
4526  /* set the next temperature */
4527  /* recompute everything at this new temp */
4528  TempChange(te[i] , true);
4529  /* converge the pressure-temperature-ionization solution for this zone */
4530  ConvPresTempEdenIoniz();
4531 
4532  /* now get the thermal emission */
4533  RT_diffuse();
4534  for(j=0;j<rfield.nflux; ++j )
4535  {
4536  agn_continuum[i][j] = rfield.ConEmitLocal[nzone][j]/(realnum)dense.eden/
4537  (dense.xIonDense[ipHYDROGEN][1] + dense.xIonDense[ipHELIUM][1] + dense.xIonDense[ipHELIUM][2] );
4538  }
4539  }
4540 
4541  /* print title for line */
4542  fprintf(ioPUN,"wl");
4543  for( i=0;i<NTE; ++i)
4544  {
4545  fprintf(ioPUN,"\tT=%.0f",te[i]);
4546  }
4547  fprintf( ioPUN , "\tcont\n");
4548 
4549  /* not print all n temperatures across a line */
4550  for(j=0;j<rfield.nflux; ++j )
4551  {
4552  fprintf( ioPUN , "%12.5e",
4553  AnuUnit(rfield.AnuOrg[j]) );
4554  /* loop over the temperatures, and for each, calculate a continuum */
4555  for( i=0;i<NTE; ++i)
4556  {
4557  fprintf(ioPUN,"\t%.3e",agn_continuum[i][j]*rfield.anu2[j]*EN1RYD/rfield.widflx[j]);
4558  }
4559  /* cont label and end of line*/
4560  fprintf( ioPUN , "\t%s\n" , rfield.chContLabel[j]);
4561  }
4562 
4563  /* now free the continua */
4564  for( i=0;i<NTE; ++i)
4565  {
4566  free( agn_continuum[i] );
4567  }
4568 
4569  /* Restore temperature stored before this routine was called */
4570  /* and force update */
4571  TempChange(TempSave , true);
4572 
4573  fprintf( ioQQQ, "AGN_Hemis - result of save AGN3 hemis - I have left the code in a disturbed state, and will now exit.\n");
4574  cdEXIT(EXIT_FAILURE);
4575 }
4576 
4577 /*SaveResults save results from save results command */
4578 /*SaveResults1Line do single line of output for the save results and save line intensity commands */
4579 STATIC void SaveResults(FILE* ioPUN)
4580 {
4581  long int i , nelem , ion;
4582 
4583  DEBUG_ENTRY( "SaveResults()" );
4584 
4585  /* used to save out line intensities, optical depths,
4586  * and column densities */
4587 
4588  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4589  input.echo(ioPUN);
4590 
4591  /* first print any cautions or warnings */
4592  cdWarnings(ioPUN);
4593  cdCautions(ioPUN);
4594  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4595 
4596  fprintf( ioPUN, "C*OPTICAL DEPTHS ELECTRON=%10.3e\n", opac.telec );
4597 
4598  fprintf( ioPUN, "BEGIN EMISSION LINES\n" );
4599  SaveResults1Line(ioPUN," ",0,0.,"Start");
4600 
4601  for( i=0; i < LineSave.nsum; i++ )
4602  {
4603  if( LineSv[i].SumLine[0] > 0. )
4604  {
4605  SaveResults1Line(ioPUN,(char*)LineSv[i].chALab,LineSv[i].wavelength,
4606  LineSv[i].SumLine[0],"Line ");
4607  }
4608  }
4609 
4610  SaveResults1Line(ioPUN," ",0,0.,"Flush");
4611 
4612  fprintf( ioPUN, " \n" );
4613 
4614  fprintf( ioPUN, "BEGIN COLUMN DENSITIES\n" );
4615 
4616  /* this dumps out the whole array,*/
4617  /* following loop relies on LIMELM being 30, assert it here in case
4618  * this is ever changed */
4619  ASSERT( LIMELM == 30 );
4620  /* this order of indices is to keep 30 as the fastest variable,
4621  * and the 32 (LIMELM+1) as the slower one */
4622  for( nelem=0; nelem<LIMELM; nelem++ )
4623  {
4624  for(ion=0; ion < nelem+1; ion++)
4625  {
4626  fprintf( ioPUN, " %10.3e", mean.xIonMean[0][nelem][ion][0] );
4627  /* throw line feed every 10 numbers */
4628  if( nelem==9|| nelem==19 || nelem==29 )
4629  {
4630  fprintf( ioPUN, "\n" );
4631  }
4632  }
4633  }
4634 
4635  fprintf( ioPUN, "END OF RESULTS\n" );
4636  fprintf( ioPUN, "**********************************************************************************************************************************\n" );
4637  return;
4638 }
4639 
4640 static const int LINEWIDTH = 6;
4641 
4642 /*SaveResults1Line do single line of output for the save results and save line intensity commands */
4643 /* the number of emission lines across one line of printout */
4644 STATIC void SaveResults1Line(
4645  FILE * ioPUN,
4646  /* 4 char + null string */
4647  const char *chLab,
4648  realnum wl,
4649  double xInten,
4650  /* null terminated string saying what to do */
4651  const char *chFunction)
4652 {
4653 
4654  long int i;
4655  static realnum wavelength[LINEWIDTH];
4656  static long ipLine;
4657  static double xIntenSave[LINEWIDTH];
4658  static char chLabSave[LINEWIDTH][5];
4659 
4660  DEBUG_ENTRY( "SaveResults1Line()" );
4661 
4662  /* if LineWidth is changed then change format in write too */
4663 
4664  if( strcmp(chFunction,"Start") == 0 )
4665  {
4666  ipLine = 0;
4667  }
4668  else if( strcmp(chFunction,"Line ") == 0 )
4669  {
4670  /* save results in array so that they can be printed when done */
4671  wavelength[ipLine] = wl;
4672  strcpy( chLabSave[ipLine], chLab );
4673  xIntenSave[ipLine] = xInten;
4674 
4675  /* now increment the counter and then check if we have filled the line,
4676  * and so should print it */
4677  ++ipLine;
4678  /* do print now if we are in column mode (one line per line) or if we have filled up
4679  * the line */
4680  if( ( strcmp(save.chPunRltType,"column") == 0 ) || ipLine == LINEWIDTH )
4681  {
4682  /* "array " is usual array 6 wide, "column" is one line per line */
4683  for( i=0; i < ipLine; i++ )
4684  {
4685  fprintf( ioPUN, " %4.4s ", chLabSave[i] );
4686  prt_wl( ioPUN, wavelength[i] );
4687  fprintf( ioPUN,"\t%.3e", xIntenSave[i]);
4688  /* >>chng 02 apr 24, do not print type */
4689  /* single column for input into data base */
4690  if( strcmp(save.chPunRltType,"column") == 0 )
4691  fprintf( ioPUN, "\n" );
4692  }
4693  /* only put cr if we did not just put one */
4694  if( strcmp(save.chPunRltType,"array ") == 0 )
4695  fprintf( ioPUN, " \n" );
4696  ipLine = 0;
4697  }
4698  }
4699  else if( strcmp(chFunction,"Flush") == 0 )
4700  {
4701  if( ipLine > 0 )
4702  {
4703  /* this is an option to print many emission lines across an output line,
4704  * the array option, or a single column of numbers, the column option
4705  * that appears on the "save results" and "save intensity" commands
4706  */
4707  /* usual array 6 wide */
4708  for( i=0; i < ipLine; i++ )
4709  {
4710  fprintf( ioPUN, " %4.4s", chLabSave[i] );
4711  prt_wl( ioPUN, wavelength[i] );
4712  fprintf( ioPUN,"\t%.3e", xIntenSave[i]);
4713  /* >>chng 02 apr 24, do not print type */
4714  /* single column for input into data base */
4715  if( strcmp(save.chPunRltType,"column") == 0 )
4716  fprintf( ioPUN, "\n" );
4717  }
4718  if( strcmp(save.chPunRltType,"array ") == 0 )
4719  fprintf( ioPUN, " \n" );
4720  }
4721  }
4722  else
4723  {
4724  fprintf( ioQQQ, " SaveResults1Line called with insane request=%5.5s\n",
4725  chFunction );
4726  cdEXIT(EXIT_FAILURE);
4727  }
4728  return;
4729 }
4730 
4731 static const int NENR_GAUNT = 37;
4732 static const int NTE_GAUNT = 21;
4733 
4734 /*SaveGaunts called by save gaunts command to output gaunt factors */
4735 STATIC void SaveGaunts(FILE* ioPUN)
4736 {
4737  long int i,
4738  charge,
4739  ite,
4740  j;
4741 
4742  realnum ener[NENR_GAUNT],
4743  ste[NTE_GAUNT],
4744  z;
4745  double tempsave;
4746  double g[NENR_GAUNT][NTE_GAUNT], gfac;
4747 
4748 
4749  DEBUG_ENTRY( "SaveGaunts()" );
4750 
4751  /* this routine is called from the PUNCH GAUNTS command
4752  * it drives the gaunt factor routine to save gaunts over full range */
4753  tempsave = phycon.te;
4754 
4755  for( i=0; i < NTE_GAUNT; i++ )
4756  {
4757  ste[i] = 0.5f*i;
4758  }
4759 
4760  for( i=0; i < NENR_GAUNT; i++ )
4761  {
4762  ener[i] = 0.5f*i - 8.f;
4763  }
4764 
4765  for( charge = 1; charge<LIMELM; charge++ )
4766  {
4767  /* nuclear charge */
4768  z = (realnum)log10((double)charge);
4769 
4770  /* energy is log of energy */
4771  for( ite=0; ite < NTE_GAUNT; ite++ )
4772  {
4773  phycon.alogte = ste[ite];
4774  phycon.te = pow(10.,phycon.alogte);
4775 
4776  for( j=0; j < NENR_GAUNT; j++ )
4777  {
4778  gfac = cont_gaunt_calc( phycon.te, (double)charge, pow(10.,(double)ener[j]) );
4779  /*fprintf(ioQQQ, "z %.2e ener %.2e temp %.2e gfac %.3e \n",
4780  pow(10.,z), pow(10.,ener[j]), (double)phycon.te, gfac );*/
4781  g[j][ite] = gfac;
4782  }
4783  }
4784 
4785  /* now save out the results */
4786  fprintf( ioPUN, " " );
4787  for( i=1; i <= NTE_GAUNT; i++ )
4788  {
4789  fprintf( ioPUN, "\t%6.1f", ste[i-1] );
4790  }
4791  fprintf( ioPUN, "\n" );
4792 
4793  for( j=0; j < NENR_GAUNT; j++ )
4794  {
4795  fprintf( ioPUN, "\t%6.2f", ener[j] );
4796  for( ite=0; ite < NTE_GAUNT; ite++ )
4797  {
4798  fprintf( ioPUN, "\t%6.2f", log10(g[j][ite]) );
4799  }
4800  fprintf( ioPUN, "\n" );
4801  }
4802 
4803  fprintf( ioPUN, " " );
4804  for( i=0; i < NTE_GAUNT; i++ )
4805  {
4806  fprintf( ioPUN, "\t%6.1f", ste[i] );
4807  }
4808  fprintf( ioPUN, "\n" );
4809 
4810  /* now do the same thing as triplets */
4811  fprintf( ioPUN, " " );
4812  for( i=0; i < NTE_GAUNT; i++ )
4813  {
4814  fprintf( ioPUN, "\t%6.1f", ste[i] );
4815  }
4816  fprintf( ioPUN, "\n" );
4817 
4818  for( i=0; i < NTE_GAUNT; i++ )
4819  {
4820  for( j=0; j < NENR_GAUNT; j++ )
4821  {
4822  fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", ste[i], ener[j],
4823  log10(g[j][i]) );
4824  }
4825  }
4826 
4827  fprintf( ioPUN, "Below is log(u), log(gamma**2), gff\n" );
4828  /* do the same thing as above, but this time print log(u) and log(gamma2) instead of temp and energy. */
4829  for( i=0; i < NTE_GAUNT; i++ )
4830  {
4831  for( j=0; j < NENR_GAUNT; j++ )
4832  {
4833  fprintf( ioPUN, "\t%6.2f\t%6.2f\t%6.2e\n", 2.*z + log10(TE1RYD) - ste[i] , log10(TE1RYD)+ener[j]-ste[i],
4834  log10(g[j][i]) );
4835  }
4836  }
4837  fprintf( ioPUN, "end of charge = %li\n", charge);
4838  fprintf( ioPUN, "****************************\n");
4839  }
4840 
4841  phycon.te = tempsave;
4842  phycon.alogte = log10( phycon.te );
4843  return;
4844 }
4845 
4846 void SaveGrid(FILE* pnunit, exit_type status)
4847 {
4848  if( pnunit == NULL )
4849  return;
4850 
4851  if( optimize.nOptimiz == 0 )
4852  {
4853  /* start of line gives abort and warning summary */
4854  fprintf( pnunit, "#Index\tFailure?\tWarnings?\tExit code\t#rank\t#seq" );
4855  /* print start of each variable command line */
4856  for( int i=0; i < grid.nintparm; i++ )
4857  {
4858  char chStr[10];
4859  strncpy( chStr, optimize.chVarFmt[i], 9 );
4860  /* make sure this small bit of string is terminated */
4861  chStr[9] = '\0';
4862  fprintf( pnunit, "\t%s", chStr );
4863  }
4864  fprintf( pnunit, "\tgrid parameter string\n" );
4865  }
4866  /* abort / warning summary for this sim */
4867  bool lgNoFailure = ( status == ES_SUCCESS || status == ES_WARNINGS );
4868  fprintf( pnunit, "%9.9ld\t%c\t%c\t%20s\t%ld\t%ld",
4869  optimize.nOptimiz,
4870  TorF(!lgNoFailure),
4871  TorF(warnings.lgWarngs),
4872  cpu.i().chExitStatus(status).c_str(),
4873  cpu.i().nRANK(),
4874  grid.seqNum );
4875  /* the grid parameters */
4876  char chGridParam[INPUT_LINE_LENGTH];
4877  char chStringHold[100];
4878  sprintf( chStringHold, "%f", grid.interpParameters[optimize.nOptimiz][0] );
4879  strcpy( chGridParam, chStringHold );
4880  for( int j=0; j < grid.nintparm; j++ )
4881  {
4882  if( j > 0 )
4883  {
4884  sprintf( chStringHold, ", %f", grid.interpParameters[optimize.nOptimiz][j] );
4885  strcat( chGridParam, chStringHold );
4886  }
4887  fprintf( pnunit, "\t%f", grid.interpParameters[optimize.nOptimiz][j] );
4888  }
4889  fprintf( pnunit, "\t%s\n", chGridParam );
4890 }
EmissionProxy::TauCon
realnum & TauCon() const
Definition: emission.h:453
chLineLbl
char * chLineLbl(const TransitionProxy &t)
Definition: transition.cpp:237
colden.h
ipFineCont
long ipFineCont(double energy_ryd)
Definition: cont_ipoint.cpp:226
GrainVar::GrainEmission
vector< realnum > GrainEmission
Definition: grainvar.h:578
thermal.h
t_pressure::PresRamCurr
double PresRamCurr
Definition: pressure.h:79
t_hyperfine::Tspin21cm
double Tspin21cm
Definition: hyperfine.h:47
t_dynamics::Heat
double Heat()
Definition: dynamics.cpp:2173
t_dense::pden
realnum pden
Definition: dense.h:98
save_line
void save_line(FILE *ip, const char *chDo, bool lgEmergent)
Definition: save_line.cpp:102
ipCOL_H2s
#define ipCOL_H2s
Definition: colden.h:18
findspecieslocal
molezone * findspecieslocal(const char buf[])
Definition: mole_species.cpp:833
TorF
char TorF(bool l)
Definition: cddefines.h:710
helike.h
lgAbort
bool lgAbort
Definition: cddefines.cpp:10
ipOXYGEN
const int ipOXYGEN
Definition: cddefines.h:312
ipT63
long ipT63
Definition: atmdat_readin.cpp:53
sprt_wl
void sprt_wl(char *chString, realnum wl)
Definition: prt.cpp:25
lines.h
h2
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
yield.h
t_phycon::EnergyBinding
double EnergyBinding
Definition: phycon.h:44
t_ionbal::PhotoRate_Shell
double **** PhotoRate_Shell
Definition: ionbal.h:111
ipHE_LIKE
const int ipHE_LIKE
Definition: iso.h:63
t_geometry::covgeo
realnum covgeo
Definition: geometry.h:35
t_dense::eden
double eden
Definition: dense.h:190
t_rfield::fine_opt_depth
realnum * fine_opt_depth
Definition: rfield.h:410
cont_gaunt_calc
double cont_gaunt_calc(double temp, double z, double photon)
Definition: cont_gaunt.cpp:26
t_dense::EdenHCorr
double EdenHCorr
Definition: dense.h:216
diatomics::H2_PunchDo
void H2_PunchDo(FILE *io, char chJOB[], const char chTime[], long int ipPun)
Definition: mole_h2_io.cpp:1224
MAX_HEADER_SIZE
static const long MAX_HEADER_SIZE
Definition: save.h:12
colliders
ColliderList colliders
Definition: collision.cpp:57
t_rfield::trans_coef_total
realnum * trans_coef_total
Definition: rfield.h:390
struc.h
dense
t_dense dense
Definition: dense.cpp:24
t_pressure::PresTurbCurr
double PresTurbCurr
Definition: pressure.h:82
lgMustPrintHeader
static bool lgMustPrintHeader
Definition: save_line.cpp:287
elementnames.h
t_optimize::nOptimiz
long int nOptimiz
Definition: optimize.h:246
t_save::lgCumulative
bool lgCumulative[LIMPUN]
Definition: save.h:219
t_mean::xIonMean
multi_arr< double, 4 > xIonMean
Definition: mean.h:14
secondaries.h
Wind::fmul
realnum fmul
Definition: wind.h:65
Singleton< t_yield >::Inst
static t_yield & Inst()
Definition: cddefines.h:175
FeIISaveLines
void FeIISaveLines(FILE *ioPUN)
Definition: atom_feii.cpp:1709
ion_recombAGN
void ion_recombAGN(FILE *io)
Definition: ion_recomb.cpp:218
t_ionbal::CotaRate
realnum CotaRate[LIMELM]
Definition: ionbal.h:242
rfield
t_rfield rfield
Definition: rfield.cpp:8
t_rfield::flux
realnum ** flux
Definition: rfield.h:86
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
geometry.h
t_pressure::pinzon
realnum pinzon
Definition: pressure.h:110
TempChange
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:51
nLevel1
long int nLevel1
Definition: taulines.cpp:28
t_hcmap::MapZone
long int MapZone
Definition: hcmap.h:20
dBaseTrans
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:17
t_pressure::pinzon_PresIntegElecThin
realnum pinzon_PresIntegElecThin
Definition: pressure.h:116
ipRecNetEsc
const int ipRecNetEsc
Definition: cddefines.h:281
diatomics::lgEnabled
bool lgEnabled
Definition: h2_priv.h:345
t_save::lg_separate_iterations
bool lg_separate_iterations[LIMPUN]
Definition: save.h:242
t_radius::Conv2PrtInten
double Conv2PrtInten
Definition: radius.h:147
t_continuum::ResolutionScaleFactor
double ResolutionScaleFactor
Definition: continuum.h:90
TexcLine
double TexcLine(const TransitionProxy &t)
Definition: transition.cpp:169
realnum
float realnum
Definition: cddefines.h:103
mole_dominant_rates
void mole_dominant_rates(const molecule *debug_species, FILE *ioOut)
Definition: mole_eval_balance.cpp:335
iterations
t_iterations iterations
Definition: iterations.cpp:5
t_elementnames::chIonStage
char chIonStage[LIMELM+1][CHARS_ION_STAGE]
Definition: elementnames.h:29
TransitionProxy::WLAng
realnum & WLAng() const
Definition: transition.h:429
t_phycon::EnergyExcitation
double EnergyExcitation
Definition: phycon.h:37
conv.h
rfield.h
t_rfield::ConSourceFcnLocal
realnum ** ConSourceFcnLocal
Definition: rfield.h:152
nWindLine
long nWindLine
Definition: cdinit.cpp:19
STATIC
#define STATIC
Definition: cddefines.h:97
t_colden::H0_21cm_lower
double H0_21cm_lower
Definition: colden.h:96
t_save::lgHashEndIter
bool lgHashEndIter[LIMPUN]
Definition: save.h:291
t_ionbal::RateRecomTot
double ** RateRecomTot
Definition: ionbal.h:198
ipExtraLymanLines
multi_arr< int, 3 > ipExtraLymanLines
Definition: taulines.cpp:24
GrainVar::TotalEden
double TotalEden
Definition: grainvar.h:528
nUTA
long int nUTA
Definition: taulines.cpp:26
t_FeII::FeIIColl
double FeIIColl[NFE2LEVN][NFE2LEVN][8]
Definition: atomfeii.h:216
t_conv::hist_temp_temp
vector< double > hist_temp_temp
Definition: conv.h:300
t_opac::TauAbsFace
realnum * TauAbsFace
Definition: opacity.h:91
ipCARBON
const int ipCARBON
Definition: cddefines.h:310
mole.h
cdColm
int cdColm(const char *chLabel, long int ion, double *theocl)
Definition: cddrive.cpp:636
molecule::index
int index
Definition: mole.h:169
t_hcmap::lgMapBeingDone
bool lgMapBeingDone
Definition: hcmap.h:33
grid
t_grid grid
Definition: grid.cpp:5
t_tag_LineSv::wavelength
realnum wavelength
Definition: lines.h:131
t_save::lgFLUSH
bool lgFLUSH
Definition: save.h:298
t_struc::testr
realnum * testr
Definition: struc.h:25
t_opac::albedo
double * albedo
Definition: opacity.h:104
t_phycon::EnergyIonization
double EnergyIonization
Definition: phycon.h:31
diatomics::H2_PunchLineStuff
void H2_PunchLineStuff(FILE *io, realnum xLimit, long index)
Definition: mole_h2_io.cpp:1167
t_dynamics::Rate
double Rate
Definition: dynamics.h:71
ipMAGNESIUM
const int ipMAGNESIUM
Definition: cddefines.h:316
t_Heavy::nsShells
long int nsShells[LIMELM][LIMELM]
Definition: heavy.h:28
t_dense::lgElmtOn
bool lgElmtOn[LIMELM]
Definition: dense.h:146
EmissionProxy::dampXvel
realnum & dampXvel() const
Definition: emission.h:553
t_save::Resolution
realnum Resolution
Definition: save.h:363
ipIRON
const int ipIRON
Definition: cddefines.h:330
t_phycon::te03
double te03
Definition: phycon.h:59
t_ionbal::RR_rate_coef_used
double ** RR_rate_coef_used
Definition: ionbal.h:212
saveFITSfile
void saveFITSfile(FILE *io, int option)
Definition: save_fits.cpp:85
t_rfield::otscon
realnum * otscon
Definition: rfield.h:195
ipHe2p3P0
const int ipHe2p3P0
Definition: iso.h:46
AnuUnit
double AnuUnit(realnum energy)
Definition: service.cpp:173
t_iso_sp::numLevels_local
long int numLevels_local
Definition: iso.h:498
t_save::ResolutionAbs
realnum ResolutionAbs
Definition: save.h:366
GrainVar::GraphiteEmission
vector< realnum > GraphiteEmission
Definition: grainvar.h:579
cpu
static t_cpu cpu
Definition: cpu.h:355
phycon
t_phycon phycon
Definition: phycon.cpp:6
t_timesc::time_H2_Dest_here
double time_H2_Dest_here
Definition: timesc.h:37
t_mole_global::list
MoleculeList list
Definition: mole.h:317
t_input::echo
void echo(FILE *ipOUT)
Definition: input.cpp:89
DoppVel
t_DoppVel DoppVel
Definition: doppvel.cpp:5
prt_LineLabels
void prt_LineLabels(FILE *ioOUT, bool lgPrintAll)
Definition: prt.cpp:168
AGN_Hemis
STATIC void AGN_Hemis(FILE *ioPUN)
Definition: save_do.cpp:4510
ipoint.h
t_save::lgSaveEveryZone
bool lgSaveEveryZone[LIMPUN]
Definition: save.h:261
ProxyIterator
Definition: proxy_iterator.h:58
HydroRecCool
double HydroRecCool(long int n, long int ipZ)
Definition: hydroreccool.cpp:10
t_rfield::outlin
realnum ** outlin
Definition: rfield.h:199
save_average
void save_average(long int ipPun)
Definition: save_average.cpp:210
t_save::wlLineList
vector< realnum > wlLineList[LIMPUN]
Definition: save.h:180
chDummy
char * chDummy
Definition: save_do.cpp:558
grid.h
ipNEON
const int ipNEON
Definition: cddefines.h:314
SDIV
sys_float SDIV(sys_float x)
Definition: cddefines.h:952
ipT370
long ipT370
Definition: atmdat_readin.cpp:42
save_colden
void save_colden(FILE *ioPUN)
Definition: save_colden.cpp:92
t_dynamics::Source
double ** Source
Definition: dynamics.h:74
t_dense::gas_phase
realnum gas_phase[LIMELM]
Definition: dense.h:71
t_radius::depth
double depth
Definition: radius.h:38
lines_service.h
EmissionProxy::Pdest
realnum & Pdest() const
Definition: emission.h:543
t_secondaries::x12tot
realnum x12tot
Definition: secondaries.h:53
t_LineSave::nsum
long int nsum
Definition: lines.h:62
TransitionProxy::ipCont
long & ipCont() const
Definition: transition.h:450
t_dense::EdenTrue
double EdenTrue
Definition: dense.h:221
t_opac::TauAbsGeo
realnum ** TauAbsGeo
Definition: opacity.h:82
Heavy
t_Heavy Heavy
Definition: heavy.cpp:5
SaveNewContinuum
STATIC void SaveNewContinuum(FILE *ioPUN)
Definition: save_do.cpp:4409
t_save::nSaveEveryZone
long int nSaveEveryZone[LIMPUN]
Definition: save.h:262
dynamics.h
t_conv::hist_pres_density
vector< double > hist_pres_density
Definition: conv.h:295
t_save::lgPunLstIter
bool lgPunLstIter[LIMPUN]
Definition: save.h:271
Wind::AccelTotalOutward
realnum AccelTotalOutward
Definition: wind.h:52
input
t_input input
Definition: input.cpp:12
ASSERT
#define ASSERT(exp)
Definition: cddefines.h:578
GrainVar::GasCoolColl
double GasCoolColl
Definition: grainvar.h:544
t_grid::interpParameters
realnum ** interpParameters
Definition: grid.h:31
t_hmi::H2_total
double H2_total
Definition: hmi.h:16
t_ionbal::UTA_ionize_rate
double ** UTA_ionize_rate
Definition: ionbal.h:170
struc
t_struc struc
Definition: struc.cpp:6
iso.h
t_colden::C1Pops
realnum C1Pops[3]
Definition: colden.h:76
t_rfield::extin_mag_V_point
double extin_mag_V_point
Definition: rfield.h:277
t_save::whichDiatomToPrint
diatomics * whichDiatomToPrint[LIMPUN]
Definition: save.h:226
t_continuum::mesh_md5sum
string mesh_md5sum
Definition: continuum.h:127
ipHYDROGEN
const int ipHYDROGEN
Definition: cddefines.h:305
t_thermal::ctot
double ctot
Definition: thermal.h:112
Wind::AccelLine
realnum AccelLine
Definition: wind.h:61
TransitionProxy::Coll
CollisionProxy Coll() const
Definition: transition.h:424
opac
t_opac opac
Definition: opacity.cpp:5
wind
Wind wind
Definition: wind.cpp:5
EXIT_SUCCESS
#define EXIT_SUCCESS
Definition: cddefines.h:138
atmdat.h
t_ionbal::DR_Badnell_rate_coef
double ** DR_Badnell_rate_coef
Definition: ionbal.h:205
t_save::nLineList
long nLineList[LIMPUN]
Definition: save.h:176
EmissionProxy::TauIn
realnum & TauIn() const
Definition: emission.h:423
NTYPES
@ NTYPES
Definition: conv.h:83
ipCOL_Hp
#define ipCOL_Hp
Definition: colden.h:26
ipoint
long ipoint(double energy_ryd)
Definition: cont_ipoint.cpp:16
t_pressure::PresIntegElecThin
realnum PresIntegElecThin
Definition: pressure.h:115
POW2
#define POW2
Definition: cddefines.h:929
GridGatherInCloudy
void GridGatherInCloudy(void)
Definition: grid_xspec.cpp:167
t_timesc::sound_speed_adiabatic
double sound_speed_adiabatic
Definition: timesc.h:45
MIN2
#define MIN2
Definition: cddefines.h:761
TransitionProxy
Definition: transition.h:23
hyperfine
t_hyperfine hyperfine
Definition: hyperfine.cpp:5
t_FeII::FeIIAul
double FeIIAul[NFE2LEVN][NFE2LEVN]
Definition: atomfeii.h:213
b1
static realnum b1[83]
Definition: atmdat_3body.cpp:25
t_phycon::te_eV
double te_eV
Definition: phycon.h:14
cddrive.h
t_species::chLabel
char * chLabel
Definition: cddefines.h:1235
DynaSave
void DynaSave(FILE *ipPnunit, char chJob)
Definition: dynamics.cpp:2148
nzone
long int nzone
Definition: cddefines.cpp:14
LineSave
t_LineSave LineSave
Definition: lines.cpp:5
FeIIPunPop
void FeIIPunPop(FILE *ioPUN, bool lgPunchRange, long int ipRangeLo, long int ipRangeHi, bool lgPunchDensity)
Definition: atom_feii.cpp:2189
FeII_Cont
realnum ** FeII_Cont
Definition: atom_feii.cpp:113
timesc
t_timesc timesc
Definition: timesc.cpp:5
cdSPEC2
void cdSPEC2(int Option, long int nEnergy, long int ipLoEnergy, long int ipHiEnergy, realnum ReturnedSpectrum[])
t_rfield::comup
double * comup
Definition: rfield.h:255
t_grid::nintparm
long nintparm
Definition: grid.h:47
radius
t_radius radius
Definition: radius.cpp:5
mole_punch
void mole_punch(FILE *punit, const char speciesname[], const char args[], bool lgHeader, bool lgData, double depth)
Definition: mole_reactions.cpp:4221
t_FeII::nFeIILevel_malloc
long int nFeIILevel_malloc
Definition: atomfeii.h:193
t_save::optname
string optname[LIMPUN]
Definition: save.h:257
t_iso_sp::nCollapsed_local
long int nCollapsed_local
Definition: iso.h:488
DynaPunchTimeDep
void DynaPunchTimeDep(FILE *ipPnunit, const char *chJob)
Definition: dynamics.cpp:2039
t_rfield::ConEmitLocal
realnum ** ConEmitLocal
Definition: rfield.h:149
t_save::emisfreq
Energy emisfreq[LIMPUN]
Definition: save.h:371
TransitionProxy::EnergyRyd
double EnergyRyd() const
Definition: transition.h:83
sexp
sys_float sexp(sys_float x)
Definition: service.cpp:914
t_save::chConPunEnr
const char * chConPunEnr[LIMPUN]
Definition: save.h:284
t_dynamics::lgTimeDependentStatic
bool lgTimeDependentStatic
Definition: dynamics.h:96
optimize
t_optimize optimize
Definition: optimize.cpp:5
t_save::lgFITS
bool lgFITS[LIMPUN]
Definition: save.h:274
t_cpu::i
t_cpu_i & i()
Definition: cpu.h:347
t_hmi::UV_Cont_rel2_Habing_TH85_depth
realnum UV_Cont_rel2_Habing_TH85_depth
Definition: hmi.h:64
TransitionProxy::Emis
EmissionList::reference Emis() const
Definition: transition.h:408
t_optimize::chVarFmt
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:263
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
t_rfield::egamry
realnum egamry
Definition: rfield.h:52
ipLine
static long int * ipLine
Definition: prt_linesum.cpp:15
t_cpu_i::big_endian
bool big_endian() const
Definition: cpu.h:287
wavelength_compare
int wavelength_compare(const void *a, const void *b)
Definition: save_do.cpp:75
t_conv::getCounterZone
long getCounterZone(const long type)
Definition: conv.h:327
dense.h
t_save::chNONSENSE
const char * chNONSENSE
Definition: save.h:233
SaveGrid
void SaveGrid(FILE *pnunit, exit_type status)
Definition: save_do.cpp:4846
ipSILICON
const int ipSILICON
Definition: cddefines.h:318
SaveLineIntensity
STATIC void SaveLineIntensity(FILE *ioPUN, long int ipPun, realnum Threshold)
Definition: save_do.cpp:4180
mole
t_mole_local mole
Definition: mole.cpp:7
TransitionProxy::Lo
qList::iterator Lo() const
Definition: transition.h:392
t_isoCTRL::SmallA
realnum SmallA
Definition: iso.h:371
ChargTranPun
void ChargTranPun(FILE *ipPnunit, char *chSave)
Definition: atmdat_char_tran.cpp:1729
t_colden::H0_21cm_upper
double H0_21cm_upper
Definition: colden.h:95
t_conv::hist_temp_cool
vector< double > hist_temp_cool
Definition: conv.h:300
FeIIPunchLevels
void FeIIPunchLevels(FILE *ioPUN)
Definition: atom_feii.cpp:1625
t_DoppVel::TurbVel
realnum TurbVel
Definition: doppvel.h:12
t_Heavy::Valence_IP_Ryd
double Valence_IP_Ryd[LIMELM][LIMELM]
Definition: heavy.h:24
cddefines.h
safe_div
sys_float safe_div(sys_float x, sys_float y, sys_float res_0by0)
Definition: cddefines.h:961
null_mole
molecule * null_mole
Definition: mole_species.cpp:64
ipSULPHUR
const int ipSULPHUR
Definition: cddefines.h:320
ipH2s
const int ipH2s
Definition: iso.h:28
t_hcmap::lgMapDone
bool lgMapDone
Definition: hcmap.h:36
FeIIPunchOpticalDepth
void FeIIPunchOpticalDepth(FILE *ioPUN)
Definition: atom_feii.cpp:1673
t_opac::opacity_sct
double * opacity_sct
Definition: opacity.h:98
t_tag_LineSv::chSumTyp
char chSumTyp
Definition: lines.h:114
t_save::lgLinEvery
bool lgLinEvery
Definition: save.h:351
ipHe1s1S
const int ipHe1s1S
Definition: iso.h:41
NTE_GAUNT
static const int NTE_GAUNT
Definition: save_do.cpp:4732
Energy::Ryd
double Ryd() const
Definition: energy.h:26
t_LineSave::ipass
long int ipass
Definition: lines.h:75
t_opac::TauScatFace
realnum * TauScatFace
Definition: opacity.h:92
t_radius::Radius
double Radius
Definition: radius.h:25
ipHe2p3P1
const int ipHe2p3P1
Definition: iso.h:47
thermal
t_thermal thermal
Definition: thermal.cpp:5
t_rfield::chLineLabel
char ** chLineLabel
Definition: rfield.h:220
ExtraLymanLines
vector< vector< TransitionList > > ExtraLymanLines
Definition: taulines.cpp:25
POW3
#define POW3
Definition: cddefines.h:936
t_hmi::H2_photodissoc_used_H2g
double H2_photodissoc_used_H2g
Definition: hmi.h:108
optimize.h
t_pressure::PresGasCurr
double PresGasCurr
Definition: pressure.h:89
TransitionProxy::Hi
qList::iterator Hi() const
Definition: transition.h:396
TotalInsanity
NORETURN void TotalInsanity(void)
Definition: service.cpp:886
t_rfield::comdn
double * comdn
Definition: rfield.h:256
SaveHeat
void SaveHeat(FILE *io)
Definition: heat_save.cpp:22
SaveResults1Line
STATIC void SaveResults1Line(FILE *ioPUN, const char *chLab, realnum wl, double xInten, const char *chFunction)
Definition: save_do.cpp:4644
plankf
double plankf(long int ip)
Definition: service.cpp:1707
lgSaveOpticalDepths
static bool lgSaveOpticalDepths
Definition: save_do.cpp:4227
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
GrainVar::bin
vector< GrainBin * > bin
Definition: grainvar.h:583
TauLine2
TransitionList TauLine2("TauLine2", &AnonStates)
t_called::lgTalk
bool lgTalk
Definition: called.h:12
ES_SUCCESS
@ ES_SUCCESS
Definition: cddefines.h:116
Wind::AccelGravity
realnum AccelGravity
Definition: wind.h:49
t_save::chHashString
char chHashString[INPUT_LINE_LENGTH]
Definition: save.h:295
nMatch
long nMatch(const char *chKey, const char *chCard)
Definition: service.cpp:451
t_tag_LineSv::SumLine
double SumLine[4]
Definition: lines.h:125
t_rfield::AnuOrg
double * AnuOrg
Definition: rfield.h:62
hyperfine.h
ipT146
long ipT146
Definition: atmdat_readin.cpp:52
Wind::dvdr
realnum dvdr
Definition: wind.h:21
t_conv::lgConvPres
bool lgConvPres
Definition: conv.h:199
t_thermal::heating
double heating[LIMELM][LIMELM]
Definition: thermal.h:158
t_pressure::PresTotlError
double PresTotlError
Definition: pressure.h:87
radius.h
t_Heavy::chShell
char chShell[7][3]
Definition: heavy.h:31
ipCOL_elec
#define ipCOL_elec
Definition: colden.h:30
FeIIPunchColden
void FeIIPunchColden(FILE *ioPUN)
Definition: atom_feii.cpp:1647
MALLOC
#define MALLOC(exp)
Definition: cddefines.h:501
colden
t_colden colden
Definition: colden.cpp:5
t_dense::xMassDensity
realnum xMassDensity
Definition: dense.h:91
t_yield::elec_eject_frac
realnum elec_eject_frac(long n, long i, long ns, long ne) const
Definition: yield.h:48
FeII
t_FeII FeII
Definition: atomfeii.cpp:5
lgFirst
static bool lgFirst
Definition: prt_linesum.cpp:16
t_geometry::lgSphere
bool lgSphere
Definition: geometry.h:24
t_rfield::nflux
long int nflux
Definition: rfield.h:43
heavy.h
t_phycon::TeProp
double TeProp
Definition: phycon.h:91
save_opacity
void save_opacity(FILE *io, long int np)
Definition: save_opacity.cpp:22
t_opac::opacity_abs
double * opacity_abs
Definition: opacity.h:95
t_save::chHeader
char chHeader[LIMPUN][MAX_HEADER_SIZE]
Definition: save.h:231
hmi.h
t_FeII::FeIINRGs
double FeIINRGs[NFE2LEVN]
Definition: atomfeii.h:207
t_iso_sp::fb
vector< freeBound > fb
Definition: iso.h:452
PrtColumns
void PrtColumns(FILE *ioMEAN, const char *chType, long int ipPun)
Definition: prt_columns.cpp:14
ipALUMINIUM
const int ipALUMINIUM
Definition: cddefines.h:317
EmissionProxy::pump
double & pump() const
Definition: emission.h:473
t_colden::colden
realnum colden[NCOLD]
Definition: colden.h:38
mean
t_mean mean
Definition: mean.cpp:17
MAX2
#define MAX2
Definition: cddefines.h:782
ionbal
t_ionbal ionbal
Definition: ionbal.cpp:5
pressure.h
t_pressure::PresTotlInit
double PresTotlInit
Definition: pressure.h:92
CollisionProxy::ColUL
realnum ColUL(const ColliderList &colls) const
Definition: collision.h:99
LIMELM
const int LIMELM
Definition: cddefines.h:258
OccupationNumberLine
double OccupationNumberLine(const TransitionProxy &t)
Definition: transition.cpp:142
RYDLAM
const UNUSED double RYDLAM
Definition: physconst.h:176
t_pressure::IntegRhoGravity
double IntegRhoGravity
Definition: pressure.h:123
cdSPEC
void cdSPEC(int Option, long int nEnergy, double ReturnedSpectrum[])
Definition: cdspec.cpp:21
ipELECTRON
@ ipELECTRON
Definition: collision.h:9
t_pressure::pres_radiation_lines_curr
double pres_radiation_lines_curr
Definition: pressure.h:101
t_conv::nPres2Ioniz
long int nPres2Ioniz
Definition: conv.h:152
t_save::lgEmergent
bool lgEmergent[LIMPUN]
Definition: save.h:216
t_hmi::HeatH2Dexc_used
double HeatH2Dexc_used
Definition: hmi.h:137
cdExecTime
double cdExecTime()
Definition: cddrive.cpp:481
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
t_iso_sp::st
qList st
Definition: iso.h:453
UTALines
TransitionList UTALines("UTALines", &AnonStates)
warnings
t_warnings warnings
Definition: warnings.cpp:11
t_rfield::getCoarseTransCoef
const realnum * getCoarseTransCoef()
Definition: rfield.cpp:10
ipRecRad
const int ipRecRad
Definition: cddefines.h:283
SaveLineData
NORETURN void SaveLineData(FILE *io)
Definition: save_linedata.cpp:25
t_save::nsave
long int nsave
Definition: save.h:222
GrainVar::dstsc
vector< double > dstsc
Definition: grainvar.h:525
SaveSpecies
void SaveSpecies(FILE *ioPUN, long int ipPun)
Definition: save_species.cpp:17
t_yield::nelec_eject
long nelec_eject(long n, long i, long ns) const
Definition: yield.h:55
t_colden::C2Pops
realnum C2Pops[5]
Definition: colden.h:64
t_isoCTRL::nLyman
long int nLyman[NISO]
Definition: iso.h:334
cdLine
long int cdLine(const char *chLabel, realnum wavelength, double *relint, double *absint)
Definition: cddrive.cpp:1228
save.h
ES_WARNINGS
@ ES_WARNINGS
Definition: cddefines.h:118
t_rfield::nupper
long int nupper
Definition: rfield.h:46
t_colden::O1Pops
realnum O1Pops[3]
Definition: colden.h:80
SaveLineStuff
STATIC void SaveLineStuff(FILE *ioPUN, const char *chJob, realnum xLimit)
Definition: save_do.cpp:4230
ipCOL_H3p
#define ipCOL_H3p
Definition: colden.h:28
t_rfield::anu2
realnum * anu2
Definition: rfield.h:79
iteration
long int iteration
Definition: cddefines.cpp:16
t_dense::xIonDense
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:125
t_save::chPunRltType
char chPunRltType[7]
Definition: save.h:319
SaveGaunts
STATIC void SaveGaunts(FILE *ioPUN)
Definition: save_do.cpp:4735
ipHe2p1P
const int ipHe2p1P
Definition: iso.h:49
t_conv::hist_pres_current
vector< double > hist_pres_current
Definition: conv.h:295
map_do
void map_do(FILE *io, const char *chType)
Definition: hcmap.cpp:23
t_save::chSave
char chSave[LIMPUN][5]
Definition: save.h:225
t_rfield::otslin
realnum * otslin
Definition: rfield.h:193
TauLines
TransitionList TauLines("TauLines", &AnonStates)
t_radius::r1r0sq
double r1r0sq
Definition: radius.h:49
t_cpu_i::lgMPI_talk
bool lgMPI_talk() const
Definition: cpu.h:328
t_save::ipPnunit
FILE * ipPnunit[LIMPUN]
Definition: save.h:197
t_save::ipConPun
long int ipConPun
Definition: save.h:287
abund
t_abund abund
Definition: abund.cpp:5
prt.h
t_save::chSaveArgs
char chSaveArgs[LIMPUN][5]
Definition: save.h:265
t_conv::hist_temp_heat
vector< double > hist_temp_heat
Definition: conv.h:300
t_save::lgRealSave
bool lgRealSave[LIMPUN]
Definition: save.h:213
t_rfield::reflin
realnum ** reflin
Definition: rfield.h:206
hydrogenic.h
EmissionProxy::Pelec_esc
realnum & Pelec_esc() const
Definition: emission.h:533
doppvel.h
t_iso_sp::RadRec_caseB
double RadRec_caseB
Definition: iso.h:513
t_rfield::fine_opac_zone
realnum * fine_opac_zone
Definition: rfield.h:408
t_grid::seqNum
long seqNum
Definition: grid.h:62
t_rfield::extin_mag_V_extended
double extin_mag_V_extended
Definition: rfield.h:281
t_thermal::htot
double htot
Definition: thermal.h:149
ipT157
long ipT157
Definition: atmdat_readin.cpp:42
t_tag_LineSv::chALab
char chALab[5]
Definition: lines.h:117
t_hmi::UV_Cont_rel2_Draine_DB96_depth
realnum UV_Cont_rel2_Draine_DB96_depth
Definition: hmi.h:74
t_rfield::nfine
long nfine
Definition: rfield.h:402
PrtMeanIon
void PrtMeanIon(char chType, bool lgDensity, FILE *)
Definition: prt_meanion.cpp:11
t_iterations::lgLastIt
bool lgLastIt
Definition: iterations.h:36
GrainVar::dstab
vector< double > dstab
Definition: grainvar.h:524
NUM_OUTPUT_TYPES
const int NUM_OUTPUT_TYPES
Definition: grid.h:21
ipH2p
const int ipH2p
Definition: iso.h:29
INPUT_LINE_LENGTH
const int INPUT_LINE_LENGTH
Definition: cddefines.h:254
grainvar.h
FeIIPunchLineStuff
void FeIIPunchLineStuff(FILE *io, realnum xLimit, long index)
Definition: atom_feii.cpp:2820
GrainVar::rate_h2_form_grains_used_total
double rate_h2_form_grains_used_total
Definition: grainvar.h:574
t_secondaries::csupra
realnum ** csupra
Definition: secondaries.h:21
FeIIPunDepart
void FeIIPunDepart(FILE *ioPUN, bool lgDoAll)
Definition: atom_feii.cpp:2029
timesc.h
rt.h
t_dense::AtomicWeight
realnum AtomicWeight[LIMELM]
Definition: dense.h:75
t_hmi::H2_Solomon_dissoc_rate_used_H2g
double H2_Solomon_dissoc_rate_used_H2g
Definition: hmi.h:92
t_conv::lgConvTemp
bool lgConvTemp
Definition: conv.h:196
t_rfield::anu
double * anu
Definition: rfield.h:58
wind.h
ionbal.h
ConvPresTempEdenIoniz
int ConvPresTempEdenIoniz(void)
Definition: conv_pres_temp_eden_ioniz.cpp:23
t_phycon::te10
double te10
Definition: phycon.h:55
magnetic.h
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
a1
static double a1[83]
Definition: atmdat_3body.cpp:27
TransitionProxy::EnergyWN
realnum & EnergyWN() const
Definition: transition.h:438
gammas.h
cdCautions
void cdCautions(FILE *ioOUT)
Definition: cddrive.cpp:226
EmissionProxy::Pesc
realnum & Pesc() const
Definition: emission.h:523
SaveDo
void SaveDo(const char *chTime)
Definition: save_do.cpp:573
t_colden::H0_ov_Tspin
double H0_ov_Tspin
Definition: colden.h:58
FindStrongestLineLabels
STATIC void FindStrongestLineLabels(void)
Definition: save_do.cpp:89
hmi
t_hmi hmi
Definition: hmi.cpp:5
is_odd
bool is_odd(int j)
Definition: cddefines.h:714
physconst.h
secondaries
t_secondaries secondaries
Definition: secondaries.cpp:5
GetDopplerWidth
realnum GetDopplerWidth(realnum massAMU)
Definition: temp_change.cpp:499
dynamics
t_dynamics dynamics
Definition: dynamics.cpp:44
t_radius::drad
double drad
Definition: radius.h:31
EmissionProxy::gf
realnum & gf() const
Definition: emission.h:513
ipCOL_H2g
#define ipCOL_H2g
Definition: colden.h:16
called
t_called called
Definition: called.cpp:5
conv
t_conv conv
Definition: conv.cpp:5
t_magnetic::pressure
double pressure
Definition: magnetic.h:33
t_tag_LineSv
Definition: lines.h:111
gv
GrainVar gv
Definition: grainvar.cpp:5
findspecies
molecule * findspecies(const char buf[])
Definition: mole_species.cpp:814
t_save::lgLineListRatio
bool lgLineListRatio[LIMPUN]
Definition: save.h:182
t_ionbal::RateIonizTot
double RateIonizTot(long nelem, long ion)
Definition: ionbal.h:254
t_save::ncSaveSkip
long int ncSaveSkip
Definition: save.h:354
PrtLinePres
void PrtLinePres(FILE *ioPRESSURE)
Definition: prt_linepres.cpp:16
EmissionProxy::damp
realnum & damp() const
Definition: emission.h:563
t_rfield::ConEmitOut
realnum ** ConEmitOut
Definition: rfield.h:161
t_save::LinEvery
long int LinEvery
Definition: save.h:350
VERSION_TRNCON
static const long VERSION_TRNCON
Definition: save.h:15
magnetic
t_magnetic magnetic
Definition: magnetic.cpp:17
t_phycon::te70
double te70
Definition: phycon.h:51
t_opac::telec
realnum telec
Definition: opacity.h:175
SaveResults
STATIC void SaveResults(FILE *ioPUN)
Definition: save_do.cpp:4579
t_pressure::PresTotlCurr
double PresTotlCurr
Definition: pressure.h:86
t_rfield::line_count
long int * line_count
Definition: rfield.h:68
dBaseSpecies
species * dBaseSpecies
Definition: taulines.cpp:14
t_rfield::widflx
realnum * widflx
Definition: rfield.h:65
Save1Line
void Save1Line(const TransitionProxy &t, FILE *ioPUN, realnum xLimit, long index, realnum DopplerWidth)
Definition: save_do.cpp:4347
t_rfield::ipV_filter
long int ipV_filter
Definition: rfield.h:259
Save_Line_RT
void Save_Line_RT(FILE *ip)
Definition: save_line.cpp:351
PrettyTransmission
realnum PrettyTransmission(long j, realnum transmission)
Definition: save_do.cpp:515
iso_ctrl
t_isoCTRL iso_ctrl
Definition: iso.cpp:6
fixit
void fixit(void)
Definition: service.cpp:991
t_rfield::ConInterOut
realnum * ConInterOut
Definition: rfield.h:164
t_FeII::FeIISTWT
double FeIISTWT[NFE2LEVN]
Definition: atomfeii.h:210
prt_wl
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:13
ipCOL_HTOT
#define ipCOL_HTOT
Definition: colden.h:12
t_rfield::outlin_noplot
realnum * outlin_noplot
Definition: rfield.h:200
t_radius::drad_x_fillfac
double drad_x_fillfac
Definition: radius.h:71
t_hmi::HeatH2Dish_used
double HeatH2Dish_used
Definition: hmi.h:129
t_phycon::alogte
double alogte
Definition: phycon.h:82
t_thermal::dCooldT
double dCooldT
Definition: thermal.h:119
ipHELIUM
const int ipHELIUM
Definition: cddefines.h:306
t_grid::lgOutputTypeOn
bool lgOutputTypeOn[NUM_OUTPUT_TYPES]
Definition: grid.h:56
t_iso_sp::trans
TransitionProxy trans(const long ipHi, const long ipLo)
Definition: iso.h:444
SaveSpecial
void SaveSpecial(FILE *io, const char *chTime)
Definition: save_special.cpp:14
taulines.h
t_phycon::telogn
double telogn[7]
Definition: phycon.h:76
RT_diffuse
void RT_diffuse(void)
Definition: rt_diffuse.cpp:34
molecule
Definition: mole.h:132
t_mole_local::species
valarray< class molezone > species
Definition: mole.h:398
t_save::chSpeciesDominantRates
char chSpeciesDominantRates[LIMPUN][CHARS_SPECIES]
Definition: save.h:368
cdWarnings
void cdWarnings(FILE *ioPNT)
Definition: cddrive.cpp:198
LINEWIDTH
static const int LINEWIDTH
Definition: save_do.cpp:4640
t_cpu_i::nRANK
long nRANK() const
Definition: cpu.h:326
t_ionbal::CollIonRate_Ground
double *** CollIonRate_Ground
Definition: ionbal.h:120
t_rfield::flux_total_incident
realnum ** flux_total_incident
Definition: rfield.h:209
t_warnings::lgWarngs
bool lgWarngs
Definition: warnings.h:56
t_radius::depth_mid_zone
double depth_mid_zone
Definition: radius.h:41
hcmap.h
ipSODIUM
const int ipSODIUM
Definition: cddefines.h:315
pressure
t_pressure pressure
Definition: pressure.cpp:5
phycon.h
geometry
t_geometry geometry
Definition: geometry.cpp:5
dBaseStates
vector< qList > dBaseStates
Definition: taulines.cpp:15
lgPopsFirstCall
static bool lgPopsFirstCall
Definition: save_do.cpp:4227
ShowMe
void ShowMe(void)
Definition: service.cpp:181
t_save::punarg
realnum punarg[LIMPUN][3]
Definition: save.h:254
EmissionProxy::TauTot
realnum & TauTot() const
Definition: emission.h:433
ipH1s
const int ipH1s
Definition: iso.h:27
t_rfield::ConEmitReflec
realnum ** ConEmitReflec
Definition: rfield.h:155
t_rfield::DiffuseLineEmission
realnum * DiffuseLineEmission
Definition: rfield.h:203
CoolSave
void CoolSave(FILE *io, char chJob[])
Definition: cool_save.cpp:20
Wind::windv
realnum windv
Definition: wind.h:18
GammaPrt
void GammaPrt(long int ipLoEnr, long int ipHiEnr, long int ipOpac, FILE *ioFILE, double total, double threshold)
Definition: cont_gammas.cpp:253
CHIANTI_Upsilon
double CHIANTI_Upsilon(long, long, long, long, double)
Definition: species2.cpp:690
iterations.h
t_DoppVel::lgTurb_pressure
bool lgTurb_pressure
Definition: doppvel.h:33
t_ADfA::ph1
realnum ph1(int i, int j, int k, int l) const
Definition: atmdat.h:329
ipT1032
long ipT1032
Definition: atmdat_readin.cpp:54
LineSvSortWL
LinSv * LineSvSortWL
Definition: cdinit.cpp:71
continuum
t_continuum continuum
Definition: continuum.cpp:5
t_phycon::EnthalpyDensity
double EnthalpyDensity
Definition: phycon.h:40
t_species::numLevels_max
long numLevels_max
Definition: cddefines.h:1239
t_rfield::emm
realnum emm
Definition: rfield.h:49
TE1RYD
const UNUSED double TE1RYD
Definition: physconst.h:183
lgCheckMonitors
bool lgCheckMonitors(FILE *ioMONITOR)
Definition: monitor_results.cpp:1603
NENR_GAUNT
static const int NENR_GAUNT
Definition: save_do.cpp:4731
AGN_He1_CS
void AGN_He1_CS(FILE *ioPun)
Definition: iso_solve.cpp:464
ipCOL_H0
#define ipCOL_H0
Definition: colden.h:22
SaveFeII_cont
STATIC realnum SaveFeII_cont(long int ipCont, long ipFeII_Cont_type)
Definition: save_do.cpp:561
opacity.h
iso_sp
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:8
t_save::chLineListLabel
vector< char * > chLineListLabel[LIMPUN]
Definition: save.h:178
called.h
mole_priv.h
t_phycon::TeInit
double TeInit
Definition: phycon.h:89
continuum.h
h2.h
LineSv
LinSv * LineSv
Definition: cdinit.cpp:70
t_save::lgPunHeader
bool lgPunHeader[LIMPUN]
Definition: save.h:246
hcmap
t_hcmap hcmap
Definition: hcmap.cpp:21
t_dynamics::lgStatic_completed
bool lgStatic_completed
Definition: dynamics.h:105
t_mole_local::findrk
double findrk(const char buf[]) const
Definition: mole_reactions.cpp:3886
EVRYD
const UNUSED double EVRYD
Definition: physconst.h:189
GrainVar::SilicateEmission
vector< realnum > SilicateEmission
Definition: grainvar.h:580
molezone::den
double den
Definition: mole.h:358
nSpecies
long int nSpecies
Definition: taulines.cpp:21
t_geometry::FillFac
realnum FillFac
Definition: geometry.h:19
nFeIIConBins
long int nFeIIConBins
Definition: atom_feii.cpp:119
t_rfield::fine_anu
realnum * fine_anu
Definition: rfield.h:412
ipHe2s3S
const int ipHe2s3S
Definition: iso.h:44
EmissionProxy::Aul
realnum & Aul() const
Definition: emission.h:613
BOLTZMANN
const UNUSED double BOLTZMANN
Definition: physconst.h:97
t_rfield::ipMaxBolt
long int ipMaxBolt
Definition: rfield.h:249
t_phycon::te
double te
Definition: phycon.h:11
t_magnetic::EnthalpyDensity
double EnthalpyDensity
Definition: magnetic.h:30
NISO
const int NISO
Definition: cddefines.h:261
monitor_results.h
ipRecEsc
const int ipRecEsc
Definition: cddefines.h:279
ipHe2p3P2
const int ipHe2p3P2
Definition: iso.h:48
warnings.h
t_opac::ipElement
long int ipElement[LIMELM][LIMELM][7][3]
Definition: opacity.h:269
t_dynamics::time_elapsed
double time_elapsed
Definition: dynamics.h:99
t_radius::Radius_mid_zone
double Radius_mid_zone
Definition: radius.h:28
t_cpu_i::chExitStatus
const string & chExitStatus(exit_type s) const
Definition: cpu.h:333
input.h
EN1RYD
const UNUSED double EN1RYD
Definition: physconst.h:179
mole_global
t_mole_global mole_global
Definition: mole.cpp:6
t_grid::lgGrid
bool lgGrid
Definition: grid.h:40
t_thermal::totcol
double totcol
Definition: thermal.h:110
atomfeii.h
Wind::AccelCont
realnum AccelCont
Definition: wind.h:55
max
long max(int a, long b)
Definition: cddefines.h:775
t_conv::hist_pres_error
vector< double > hist_pres_error
Definition: conv.h:295
t_mole_global::num_calc
int num_calc
Definition: mole.h:314
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
ipH_LIKE
const int ipH_LIKE
Definition: iso.h:62
save
t_save save
Definition: save.cpp:5
t_rfield::TotDiff2Pht
realnum * TotDiff2Pht
Definition: rfield.h:187
t_rfield::ConRefIncid
realnum ** ConRefIncid
Definition: rfield.h:167
t_dynamics::Cool
double Cool()
Definition: dynamics.cpp:2187
t_save::ipEmisFreq
long ipEmisFreq[LIMPUN]
Definition: save.h:372
t_pressure::PresInteg
realnum PresInteg
Definition: pressure.h:109
HFLines
TransitionList HFLines("HFLines", &AnonStates)
SMALLFLOAT
const realnum SMALLFLOAT
Definition: cpu.h:191
exit_type
exit_type
Definition: cddefines.h:115
ipT610
long ipT610
Definition: atmdat_readin.cpp:41
wavelength
static realnum * wavelength
Definition: monitor_results.cpp:70
TransitionList::Emis
EmissionList & Emis()
Definition: transition.h:329
t_rfield::chContLabel
char ** chContLabel
Definition: rfield.h:223
mean.h
g
static double * g
Definition: species2.cpp:28
Generated by   doxygen 1.8.17