cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
save_linedata.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2017 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*SaveLineData punches selected line data for all lines transferred in code */
4 /*Save1LineData save data for one line */
5 #include "cddefines.h"
6 #include "taulines.h"
7 #include "iso.h"
8 #include "phycon.h"
9 #include "elementnames.h"
10 #include "dense.h"
11 #include "conv.h"
12 #include "lines.h"
13 #include "prt.h"
14 #include "h2.h"
15 #include "thermal.h"
16 #include "cooling.h"
17 #include "save.h"
18 #include "mole.h"
19 
20 NORETURN void SaveLineData(FILE * ioPUN)
21 {
22  const long nskip=2; /* number of emission lines per line of output */
23  double tot;
24  bool lgElemOff=false;
25 
26  DEBUG_ENTRY( "SaveLineData()" );
27 
28  /* routine punches out (on unit ioPUN) line data
29  * for all recombination lines, and all transitions that are transferred */
30 
31  /* say what is happening so we know why we stopped */
32  fprintf( ioQQQ, " saving line data, then stopping\n" );
33 
34  /* first check that all lines are turned on */
35  for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
36  {
37  if( !dense.lgElmtOn[nelem] )
38  {
39  fprintf(ioQQQ," WARNING - I am saving line data but element %s is turned off.\n",
41  lgElemOff = true;
42  }
43  }
44  if( lgElemOff )
45  {
46  fprintf(ioQQQ,"Some elements are turned off and save line data requested.\n");
47  fprintf(ioQQQ,"Code is now designed to do save line data only with all elements on.\n");
48  fprintf(ioQQQ,"Please try again with all elements on.\n");
49  fprintf(ioQQQ,"Please try again with all elements on.\n");
51  }
52 
53  /* evaluate rec coefficient at constant temperature if this is set, else
54  * use 10,000K */
55  double TeNew;
57  {
58  TeNew = thermal.ConstTemp;
59  }
60  else
61  {
62  TeNew = 1e4;
63  }
64  TempChange(TeNew , false);
65 
66  /* this is set of Dima's recombination lines */
68  fprintf( ioPUN, "\n Recombination lines of C, N, O\n" );
69  fprintf( ioPUN, "#Ion\tWL(A)\tCoef\tIon\tWL(A)\tCoef\n" );
70  for( long i=0; i<NRECCOEFCNO; i+=nskip)
71  {
72  /* nskip is set to 2 above */
73  long limit = MIN2(NRECCOEFCNO,i+nskip);
74  fprintf( ioPUN, " " );
75  for( long j=i; j < limit; j++ )
76  {
77  string chLabel = chIonLbl( LineSave.RecCoefCNO[0][i], long(LineSave.RecCoefCNO[0][i]-LineSave.RecCoefCNO[1][i]+1.01) );
78 
79  fprintf( ioPUN, "%s\t%6ld\t%8.3f\t",
80  chLabel.c_str(),
81  (long)(LineSave.RecCoefCNO[2][j]+0.5),
82  log10(SDIV(LineSave.RecCoefCNO[3][j]) ) );
83  }
84  fprintf( ioPUN, " \n" );
85  }
86  fprintf( ioPUN, "\n\n" );
87 
89  dense.EdenHCorr = 1.;
90  EdenChange( 1. );
91 
92  /* want very small neutral fractions so get mostly e- cs */
93  dense.xIonDense[ipHYDROGEN][1] = 1.e-5f;
94  findspecieslocal("H2")->den = 0.;
95  dense.xIonDense[ipHYDROGEN][1] = 1.;
96 
97  for( long i=0; i < nWindLine; i++ )
98  {
99  TauLine2[i].Lo()->Pop() = 1.;
100  }
101 
102  for( size_t i=0; i < UTALines.size(); i++ )
103  {
104  (*UTALines[i].Lo()).Pop() = 1.;
105  }
106 
107  for( long i=0; i < LIMELM; i++ )
108  {
109  for( long j=0; j < LIMELM+1; j++ )
110  {
111  dense.xIonDense[i][j] = 1.;
112  }
113  }
114 
115  /* evaluate cooling, this forces evaluation of collision strengths */
116  CoolEvaluate(&tot);
117 
118  fprintf( ioPUN, " Level 2 transferred lines\n" );
119  PrintLineDataHeader( ioPUN );
120  for( long i=0; i < nWindLine; i++ )
121  {
122  if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
123  {
124  Save1LineData( TauLine2[i] , ioPUN , true );
125  }
126  }
127 
128  fprintf( ioPUN, "\n\n\n end level 2, start inner shell UTA\n" );
129  PrintLineDataHeader( ioPUN );
130  for( size_t i=0; i < UTALines.size(); i++ )
131  {
132  Save1LineData( UTALines[i] , ioPUN , true );
133  }
134 
135  fprintf( ioPUN, "\n\n\n end inner shell, start H-like iso seq\n" );
136  /* h-like iso sequence */
137  /* the hydrogen like iso-sequence */
138  PrintLineDataHeader( ioPUN );
139  for( long nelem=0; nelem < LIMELM; nelem++ )
140  {
141  if( dense.lgElmtOn[nelem] )
142  {
143  iso_collide( ipH_LIKE, nelem );
144  /* arrays are dim'd iso_sp[ipH_LIKE][nelem].numLevels_max+1 */
145  /* keep this limit to iso.numLevels_max, instead of _local. */
146  for( long ipLo=ipH1s; ipLo < iso_sp[ipH_LIKE][nelem].numLevels_max-1; ipLo++ )
147  {
148  for( long ipHi=ipLo+1; ipHi < iso_sp[ipH_LIKE][nelem].numLevels_max; ipHi++ )
149  {
150  Save1LineData( iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo) , ioPUN , false );
151  }
152  }
153  }
154  }
155 
156  fprintf( ioPUN, "\n\n\n end H-like iso seq, start He-like iso seq\n" );
157  PrintLineDataHeader( ioPUN );
158  for( long nelem=1; nelem < LIMELM; nelem++ )
159  {
160  if( nelem < 2 || dense.lgElmtOn[nelem] )
161  {
162  /* arrays are dim'd iso_sp[ipH_LIKE][nelem].numLevels_max+1 */
163  for( long ipLo=ipHe1s1S; ipLo < iso_sp[ipHE_LIKE][nelem].numLevels_max-1; ipLo++ )
164  {
165  for( long ipHi=ipLo+1; ipHi < iso_sp[ipHE_LIKE][nelem].numLevels_max; ipHi++ )
166  {
167  Save1LineData( iso_sp[ipHE_LIKE][nelem].trans(ipHi,ipLo) , ioPUN , false );
168  }
169  }
170  }
171  }
172 
173  fprintf( ioPUN, "\n\n\n end he-like iso seq, start hyperfine structure lines\n" );
174  /* fine structure lines */
175  PrintLineDataHeader( ioPUN );
176  for( size_t i=0; i < HFLines.size(); i++ )
177  {
178  Save1LineData( HFLines[i] , ioPUN , true );
179  }
180 
181  fprintf( ioPUN, "\n\n\n end hyperfine, start database lines\n" );
182  /* Databases: Atoms & Molecules*/
183  PrintLineDataHeader( ioPUN );
184  for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
185  {
186  for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
187  em != dBaseTrans[ipSpecies].Emis().end(); ++em)
188  {
189  Save1LineData( (*em).Tran() , ioPUN , true );
190  }
191  }
192 
193  fprintf( ioPUN, "\n\n\n end database, start satellite lines\n" );
194  PrintLineDataHeader( ioPUN );
195  for( long ipISO = ipHE_LIKE; ipISO < NISO; ipISO++ )
196  {
197  for( long nelem = ipISO; nelem < LIMELM; nelem++ )
198  {
199  if( dense.lgElmtOn[nelem] && iso_ctrl.lgDielRecom[ipISO] )
200  {
201  for( long i=0; i<iso_sp[ipISO][nelem].numLevels_max; i++ )
202  {
203  Save1LineData( SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][i]],
204  ioPUN , true );
205  }
206  }
207  }
208  }
209 
210  /* want very small ionized fractions so get mostly H2 cs */
212  dense.EdenHCorr = 1e-6f;
213  findspecieslocal("H2")->den = 1.;
214  findspecieslocal("H2*")->den = 1.;
215  dense.xIonDense[ipHYDROGEN][1] = 1e-6;
216  EdenChange( 1e-6 );
217 
218  /* H2 molecule */
219  fprintf( ioPUN, "\n\n\n" );
220  fprintf( ioPUN, " end satellite, start H2 lines\n" );
221 
222  /* ioPUN unit, and option to print all possible lines - false indicates
223  * save only significant lines */
224  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
225  {
226  bool lgPopsConverged;
227  double old_val, new_val;
228  (*diatom)->H2_LevelPops( lgPopsConverged, old_val, new_val );
229  (*diatom)->H2_Punch_line_data( ioPUN, false );
230  }
231 
232  fprintf( ioPUN, "\n\n\n" );
233  fprintf( ioPUN, " end H2\n" );
234 
235  /* ChkMonitorstring is searched for by one of the scripts in the nightly run
236  * this run will terminate with no asserts but that is the correct behavior */
237  fprintf( ioQQQ , "\n The code is left in a disturbed state after creating the SAVE LINE DATA file.\n"
238  " No calculation is actually performed, only the SAVE LINE DATA file is produced.\n"
239  " Remove the SAVE LINE DATA command to do the calculation.\n\n ChkMonitorend is ok.\n" );
240 
241  /* stop when done, we have done some serious damage to the code */
243 }
244 
245 /* Print header for Save1LineData() function */
246 void PrintLineDataHeader( FILE * ioPUN )
247 {
248  fprintf( ioPUN, "#Ion\tWL\tgl\tgu\tgf\tA\tCS\tn(crt)\tdamp\n" );
249 }
250 
251 
252 /*Save1LineData save data for one line */
254  const TransitionProxy& t,
255  FILE * ioPUN,
256  /* flag saying whether to give collision strength too - in multi level atoms
257  * it will be not valid without a great deal more work */
258  bool lgCS_2 )
259 {
260  double CritDen;
261 
262  DEBUG_ENTRY( "Save1LineData()" );
263 
264  if( t.ipCont() <= 0 )
265  {
266  // not a real line, just give \n
267  return;
268  }
269 
275  /*iWL = iWavLen( t , &chUnits , &chShift );*/
276  /* ion label, like C 1 */
277  fprintf(ioPUN,"%s\t", chIonLbl( t ).c_str() );
278 
279  /* this is the second piece of the line label, pick up string after start */
280 
281  /* the wavelength */
282  if( strcmp( save.chConSavEnr[save.ipConPun], "labl" )== 0 )
283  {
284  prt_wl( ioPUN , t.WLAng() );
285  }
286  else
287  {
288  /* this converts energy in Rydbergs into any of the other units */
289  fprintf( ioPUN , "%.5e", AnuUnit((realnum)(t.EnergyRyd())) );
290  }
291 
292  fprintf( ioPUN, "\t%3ld\t%3ld",
293  /* lower and upper stat weights */
294  (long)((*t.Lo()).g()),
295  (long)((*t.Hi()).g()) );
296 
297  /* oscillator strength */
298  fprintf( ioPUN,PrintEfmt("\t%9.2e", t.Emis().gf()));
299 
300  /* Einstein A for transition */
301  fprintf( ioPUN,PrintEfmt("\t%9.2e", t.Emis().Aul()));
302 
303  /* next collision strengths, use different formats depending on size
304  * of collision strength */
305  if( t.Coll().col_str() > 100. )
306  {
307  fprintf( ioPUN, "\t%7.1f", t.Coll().col_str() );
308  }
309  else if( t.Coll().col_str() > 10. )
310  {
311  fprintf( ioPUN, "\t%7.2f", t.Coll().col_str() );
312  }
313  else if( t.Coll().col_str() > 1. )
314  {
315  fprintf( ioPUN, "\t%7.3f", t.Coll().col_str() );
316  }
317  else if( t.Coll().col_str() > .01 )
318  {
319  fprintf( ioPUN, "\t%7.4f", t.Coll().col_str() );
320  }
321  else if( t.Coll().col_str() > 0.0 )
322  {
323  fprintf( ioPUN, "\t%.3e", t.Coll().col_str() );
324  }
325  else
326  {
327  fprintf( ioPUN, "\t%7.4f", 0. );
328  }
329 
330  /* now print critical density but only if cs is positive
331  * >>chng 06 mar 24, add flag lgCS_2 - in multi-level systems do not want
332  * to save cs since not computed properly */
333  if( lgCS_2 && t.Coll().col_str()> 0. )
334  {
335  CritDen = t.Emis().Aul() * (*t.Hi()).g()*phycon.sqrte / (t.Coll().col_str()*COLL_CONST);
336  }
337  else
338  {
339  CritDen = 0.;
340  }
341  fprintf( ioPUN, "\t%.3e",CritDen );
342 
343  // damping constant for current conditions
344  fprintf( ioPUN,PrintEfmt("\t%9.2e", t.Emis().damp() ));
345 
346  fprintf( ioPUN, "\n" );
347 
348  return;
349 }
void prt_wl(FILE *ioOUT, realnum wl)
Definition: prt.cpp:44
t_thermal thermal
Definition: thermal.cpp:6
void SetGasPhaseDensity(const long nelem, const realnum density)
Definition: dense.cpp:106
string chIonLbl(const TransitionProxy &t)
Definition: transition.cpp:230
size_t size(void) const
Definition: transition.h:331
TransitionList UTALines("UTALines",&AnonStates)
const int ipHE_LIKE
Definition: iso.h:65
#define NORETURN
Definition: cpu.h:451
double EdenHCorr
Definition: dense.h:227
multi_arr< int, 3 > ipSatelliteLines
Definition: taulines.cpp:34
t_isoCTRL iso_ctrl
Definition: iso.cpp:9
const int NISO
Definition: cddefines.h:310
void Save1LineData(const TransitionProxy &t, FILE *io, bool lgCS_2)
const int NRECCOEFCNO
Definition: atmdat_adfa.h:7
#define PrintEfmt(F, V)
Definition: cddefines.h:1364
TransitionList HFLines("HFLines",&AnonStates)
t_phycon phycon
Definition: phycon.cpp:6
t_LineSave LineSave
Definition: lines.cpp:9
void CoolEvaluate(double *tot)
Definition: cool_eval.cpp:58
NORETURN void SaveLineData(FILE *io)
bool lgTemperatureConstant
Definition: thermal.h:44
FILE * ioQQQ
Definition: cddefines.cpp:7
molezone * findspecieslocal(const char buf[])
const int ipHe1s1S
Definition: iso.h:43
TransitionList TauLine2("TauLine2",&AnonStates)
#define MIN2(a, b)
Definition: cddefines.h:807
long int nSpecies
Definition: taulines.cpp:22
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:31
t_dense dense
Definition: global.cpp:15
static t_ADfA & Inst()
Definition: cddefines.h:209
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
realnum & gf() const
Definition: emission.h:558
const int ipH1s
Definition: iso.h:29
EmissionList::reference Emis() const
Definition: transition.h:447
long & ipCont() const
Definition: transition.h:489
bool lgDielRecom[NISO]
Definition: iso.h:385
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
vector< diatomics * > diatoms
Definition: h2.cpp:8
qList::iterator Hi() const
Definition: transition.h:435
#define cdEXIT(FAIL)
Definition: cddefines.h:484
vector< vector< TransitionList > > SatelliteLines
Definition: taulines.cpp:35
long nWindLine
Definition: cdinit.cpp:19
double AnuUnit(realnum energy)
Definition: service.cpp:197
bool lgElmtOn[LIMELM]
Definition: dense.h:160
qList::iterator Lo() const
Definition: transition.h:431
realnum ConstTemp
Definition: thermal.h:56
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:307
double den
Definition: mole.h:421
realnum & col_str() const
Definition: collision.h:191
realnum RecCoefCNO[4][NRECCOEFCNO]
Definition: lines.h:126
CollisionProxy Coll() const
Definition: transition.h:463
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double EnergyRyd() const
Definition: transition.h:95
realnum & damp() const
Definition: emission.h:608
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1217
const char * chConSavEnr[LIMPUN]
Definition: save.h:387
sys_float SDIV(sys_float x)
Definition: cddefines.h:1006
long int numLevels_max
Definition: iso.h:524
void rec_lines(double t, realnum r[][NRECCOEFCNO])
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
void EdenChange(double EdenNew)
Definition: eden_change.cpp:11
double sqrte
Definition: phycon.h:58
vector< TransitionList > dBaseTrans
Definition: taulines.cpp:18
t_save save
Definition: save.cpp:5
double te
Definition: phycon.h:21
const int ipHYDROGEN
Definition: cddefines.h:348
realnum & Aul() const
Definition: emission.h:668
void iso_collide(long ipISO, long nelem)
long int ipConPun
Definition: save.h:390
realnum & WLAng() const
Definition: transition.h:468
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
void PrintLineDataHeader(FILE *ioPUN)
#define EXIT_SUCCESS
Definition: cddefines.h:166