cloudy  trunk
parse_atom_h2.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 /*ParseAtomH2 parse information from the atom command line */
4 #include "cddefines.h"
5 #include "hmi.h"
6 #include "h2.h"
7 #include "h2_priv.h"
8 #include "parser.h"
9 #include "thirdparty.h"
10 #include "deuterium.h"
11 
12 /*ParseAtomH2 parse information from the atom command line */
13 void ParseAtomH2(Parser &p )
14 {
15  long int j;
16 
17  DEBUG_ENTRY( "ParseAtomH2()" );
18 
19  fixit(); // this must be generalized!!!
20  // easiest way is to create and populate diatoms before parsing,
21  // so we can simply iterator of diatoms and match strings.
22  // probably will want to put label on command line in quotes
23  // and change command to, for example,
24  // diatom "H2"
25  // do that work before entering this routine.
26  diatomics *diatom = NULL;
27 
28  if( p.nMatch(" H2 " ) )
29  {
30  diatom = &h2;
31  /* this command has a 2 in the H2 label - must not parse the two by
32  * accident. Get the first number off the line image, and confirm that
33  * it is a 2 */
34  j = (long int)p.FFmtRead();
35  if( j != 2 )
36  {
37  fprintf( ioQQQ, " Something is wrong with the order of the numbers on this line.\n" );
38  fprintf( ioQQQ, " The first number I encounter should be a 2.\n Sorry.\n" );
39  cdEXIT(EXIT_FAILURE);
40  }
41  }
42 
43  /* the mere calling of this routine turns the large H2 molecule on */
44  diatom->lgEnabled = true;
45 
46  if( p.nMatch("LEVE") )
47  {
48  /* number of electronic levels */
49 
50  /* lgREAD_DATA is false at start of calculation, set true when
51  * space allocated for the H lines. Once done we must ignore all
52  * future changes in the number of levels */
53  if( !diatom->lgREAD_DATA )
54  {
55  diatom->n_elec_states = (long int)p.FFmtRead();
56  if( p.lgEOL() )
57  {
58  if( p.nMatch("LARG") )
59  {
60  /* LARGE is option to use the most number of electronic levels */
61  diatom->n_elec_states = N_ELEC;
62  }
63  else
64  {
65  p.NoNumb("number of electronic levels");
66  }
67  }
68 
69  /* do not allow fewer than 3 - that includes Lyman & Werner bands */
70  if( diatom->n_elec_states < 3 )
71  {
72  fprintf( ioQQQ, " This would be too few electronic levels - resetting to 3.\n" );
73  diatom->n_elec_states = 3;
74  }
75  /* N_ELEC is the greatest number of elec lev possible */
76  else if( diatom->n_elec_states > N_ELEC )
77  {
78  fprintf( ioQQQ,
79  " This would be too many levels, the limit is %i.\n" ,
80  N_ELEC);
81  cdEXIT(EXIT_FAILURE);
82  }
83  }
84  }
85 
86  else if( p.nMatch("LIMI") )
87  {
88  /* the limit to the H2 / Htot ratio -
89  * if smaller than this, do not compute large H2 mole */
90  diatom->H2_to_H_limit = p.FFmtRead();
91  if( p.lgEOL() )
92  {
93  /* did not find a number, either mistake or key " off" */
94  if( p.nMatch( " OFF" ) )
95  {
96  /* turn off limit */
97  diatom->H2_to_H_limit = -1.;
98  }
99  else
100  {
101  p.NoNumb( "limit to the H2 / Htot ratio" );
102  }
103  }
104  else
105  {
106  /* got a number, check if negative and so a log */
107  /* a number <= 0 is the log of the ratio */
108  if( diatom->H2_to_H_limit <= 0. )
109  diatom->H2_to_H_limit = pow(10., diatom->H2_to_H_limit);
110  }
111  }
112  else if( p.nMatch("GBAR" ) )
113  {
114  /* option to either use, or not use, gbar approximation for low X
115  * levels with no collision data - by default it is on */
116  if( p.nMatch(" OFF" ) )
117  {
118  diatom->lgColl_gbar = false;
119  }
120  else if( p.nMatch(" ON " ) )
121  {
122  diatom->lgColl_gbar = true;
123  }
124  else
125  {
126  fprintf( ioQQQ,
127  " The gbar approximation must be off (\" OFF\") or on (\" ON \").\n");
128  cdEXIT(EXIT_FAILURE);
129  }
130  }
131  /* option to turn collisional effects off or on */
132  else if( p.nMatch("COLL" ) )
133  {
134  /* option to turn collisional dissociation off or on */
135  if( p.nMatch("DISS" ) )
136  {
137  /* option to turn collisions off */
138  if( p.nMatch(" ON " ) )
139  {
140  /* this is the default, leave collisions off */
141  diatom->lgColl_dissoc_coll = true;
142  }
143  else
144  {
145  /* default (and only reason for this command) is to turn off collisions */
146  diatom->lgColl_dissoc_coll = false;
147  }
148  }
149  /* option to turn collisional dissociation off or on
150  * >>chng 06 mar 01, had been simply if - so all collisions were turned off
151  * when dissociation collisions turned off -
152  * due to bucket else at end */
153  else if( p.nMatch("ORTH" ) && p.nMatch("PARA" ) )
154  {
155  /* option to turn ortho - para collisions with particles off */
156  if( p.nMatch(" ON " ) )
157  {
158  /* this is the default, leave collisions off */
159  diatom->lgH2_ortho_para_coll_on = true;
160  }
161  else
162  {
163  /* default (and only reason for this command) is to turn off
164  * ortho-para collisions */
165  diatom->lgH2_ortho_para_coll_on = false;
166  }
167  }
168 
169  /* option to turn collisional effects off or on */
170  else if( p.nMatch("GRAI" ) )
171  {
172  /* option to turn collisions off */
173  if( p.nMatch(" ON" ) )
174  {
175  /* this is the default, leave collisions off */
176  diatom->lgH2_grain_deexcitation = true;
177  }
178  else
179  {
180  /* default (and only reason for this command) is to turn off collisions */
181  diatom->lgH2_grain_deexcitation = false;
182  }
183  }
184  else if( p.nMatch(" HE " ) )
185  {
186  /* atom H2 He collisions ORNL (the default), Le BOURlot, and OFF
187  * which data set for He collisions,
188  * Teck Lee et al. ApJ to be submitted */
189  if( p.nMatch(" NEW" ) || p.nMatch("ORNL" ) )
190  {
191  /* use the new coefficients */
192  diatom->lgH2_He_ORNL = true;
193  diatom->coll_source[1].filename = "coll_rates_He_ORNL.dat";
194  }
195  else if( p.nMatch(" OLD" ) || p.nMatch("BOUR" ) )
196  {
197  /* use the coefficients from
198  *>>refer H2 collision Le Bourlot, J., Pineau des Forets,
199  *>>refercon G., & Flower, D.R. 1999, MNRAS, 305, 802*/
200  diatom->lgH2_He_ORNL = false;
201  diatom->coll_source[1].filename = "coll_rates_He_LeBourlot.dat";
202  }
203  else
204  {
205  fprintf( ioQQQ,
206  " I did not find a keyword on this ATOM H2 HE command - I know about the keys ORNL and Le BOURlot\n");
207  cdEXIT(EXIT_FAILURE);
208  }
209  }
210 
211  /*>>chng 08 feb 27, GS*/
212  else if( p.nMatch("ORH2" ) )
213  {
214  /* atom H2 H2ortho collisions ORNL (the default), Le BOURlot, and OFF
215  * which data set for H2 collisions,
216  * Teck Lee et al. ApJ to be submitted */
217  if( p.nMatch("ORNL" ) )
218  {
219  /* use the new coefficients */
220  diatom->lgH2_ORH2_ORNL = true;
221  diatom->coll_source[2].filename = "coll_rates_H2ortho_ORNL.dat";
222  }
223  else if( p.nMatch("BOUR" ) )
224  {
225  /* use the coefficients from
226  *>>refer H2 collision Le Bourlot, J., Pineau des Forets,
227  *>>refercon G., & Flower, D.R. 1999, MNRAS, 305, 802*/
228  diatom->lgH2_ORH2_ORNL = false;
229  diatom->coll_source[2].filename = "coll_rates_H2ortho_LeBourlot.dat";
230  }
231  else
232  {
233  fprintf( ioQQQ,
234  " I did not find a keyword on this ATOM H2 ohH2 command - I know about the keys ORNL and Le BOURlot\n");
235  cdEXIT(EXIT_FAILURE);
236  }
237  }
238 
239  else if( p.nMatch("PAH2" ) )
240  {
241  /* atom H2 H2ortho collisions ORNL (the default), Le BOURlot, and OFF
242  * which data set for He collisions,
243  * Teck Lee et al. ApJ to be submitted */
244  if( p.nMatch("ORNL" ) )
245  {
246  /* use the new coefficients */
247  diatom->lgH2_PAH2_ORNL = true;
248  diatom->coll_source[3].filename = "coll_rates_H2para_ORNL.dat";
249  }
250  else if( p.nMatch("BOUR" ) )
251  {
252  /* use the coefficients from
253  *>>refer H2 collision Le Bourlot, J., Pineau des Forets,
254  *>>refercon G., & Flower, D.R. 1999, MNRAS, 305, 802*/
255  diatom->lgH2_PAH2_ORNL = false;
256  diatom->coll_source[3].filename = "coll_rates_H2para_LeBourlot.dat";
257  }
258  else
259  {
260  fprintf( ioQQQ,
261  " I did not find a keyword on this ATOM H2 paH2 command - I know about the keys ORNL and Le BOURlot\n");
262  cdEXIT(EXIT_FAILURE);
263  }
264  }
265 
266  else
267  {
268  /* option to turn all collisions off */
269  if( p.nMatch(" ON " ) )
270  {
271  /* this is the default, leave collisions on */
272  diatom->lgColl_deexec_Calc = true;
273  }
274  else
275  {
276  /* default (and only reason for this command) is to turn off collisions */
277  diatom->lgColl_deexec_Calc = false;
278  }
279  }
280  }
281 
282  /* set number of levels in matrix, but not trace matrix option */
283  else if( p.nMatch("MATR" ) && !p.nMatch("TRAC" ) )
284  {
285  /* matrix option sets the number of levels that will
286  * be included in the matrix solution */
287  long numLevels = (long)p.FFmtRead();
288  if( p.nMatch(" ALL") )
289  {
290  /* " all" means do all of X, but space has not yet been allocated,
291  * so we do know know how many levels are within X - set special
292  * flag that will be used then this is known */
293  numLevels = -1;
294  }
295  else if( p.lgEOL() && !(p.nMatch(" OFF") || p.nMatch("NONE") ) )
296  {
297  /* this branch hit eol but OFF or NONE is not on line - this is a mistake */
298  fprintf( ioQQQ,
299  " The total number of levels used in the matrix solver must be entered, or keywords ALL or NONE entered.\n Sorry.\n");
300  cdEXIT(EXIT_FAILURE);
301  }
302 
303  diatom->set_numLevelsMatrix( numLevels );
304  /* cannot check less than total number of levels within x since not yet set
305  * We do not certify that matrix limits are greater than 1 -
306  * zero or <0 limits just turns if off, as did the off option */
307  }
308  else if( p.nMatch(" LTE" ) )
309  {
310  /* LTE option causes code to assume LTE for level populations */
311  diatom->lgLTE = true;
312  }
313 
314  else if( p.nMatch("TRAC" ) )
315  {
316  /* these are used to set trace levels of output
317  diatom->n_trace_final = 1;
318  diatom->n_trace_iterations = 2;
319  diatom->n_trace_full = 3;
320  diatom->n_trace_matrix = 4*/
321 
322  /* turns on trace printout - there are multiple levels */
323  if( p.nMatch("FINA" ) )
324  {
325  /* FINAL gives only final information when solver exits */
326  diatom->nTRACE = diatom->n_trace_final;
327  }
328  else if( p.nMatch("ITER" ) )
329  {
330  /* follow iterations within each call */
331  diatom->nTRACE = diatom->n_trace_iterations;
332  }
333  else if( p.nMatch("FULL" ) )
334  {
335  /* full details of solution - this is also the default*/
336  diatom->nTRACE = diatom->n_trace_full;
337  }
338  else if( p.nMatch("MATR" ) )
339  {
340  /* print the matrices used for X */
341  diatom->nTRACE = diatom->n_trace_matrix;
342  }
343  else
344  {
345  /* full details of solution is also the default*/
346  diatom->nTRACE = diatom->n_trace_full;
347  }
348  }
349  else if( p.nMatch("NOIS" ) )
350  {
351  unsigned int iseed;
352  /* check on effects of uncertainties in collision rates */
353  diatom->lgH2_NOISE = true;
354  diatom->lgH2_NOISECOSMIC = true;
355 
356  /* optional mean - default is 0 */
357  diatom->xMeanNoise = p.FFmtRead();
358  if( p.lgEOL() )
359  diatom->xMeanNoise = 0.;
360 
361  /* this is the standard deviation for the mole, with default */
362  diatom->xSTDNoise = p.FFmtRead();
363  if( p.lgEOL() )
364  diatom->xSTDNoise = 0.5;
365 
366  /* this may be a seed for the random number generator. if no seed is
367  * set then use system time, and always get different sequence */
368  iseed = (unsigned int)p.FFmtRead();
369  /* returned 0 if eol hit */
370  if( iseed > 0 )
371  {
372  /* user set seed */
373  init_genrand( iseed );
374  }
375  else
376  {
377  init_genrand( (unsigned)time( NULL ) );
378  }
379  }
380 
381  else if( p.nMatch("THER" ) )
382  {
383  /* change the treatment of the heating - cooling effects of H2,
384  * options are simple (use TH85 expressions) and full (use large molecule)*/
385  if( p.nMatch("SIMP" ) )
386  {
387  hmi.lgH2_Thermal_BigH2 = false;
388  }
389  else if( p.nMatch("FULL" ) )
390  {
391  /* this is the default - use big atom */
392  hmi.lgH2_Thermal_BigH2 = true;
393  }
394  }
395 
396  else if( p.nMatch("CHEM" ) )
397  {
398  /* atom h2 chemistry simple command
399  * change the treatment of the chemistry - formation and destruction,
400  * options are simple (use TH85 expressions) and full (use large molecule)*/
401  if( p.nMatch("SIMP" ) )
402  {
403  hmi.lgH2_Chemistry_BigH2 = false;
404  }
405  else if( p.nMatch("FULL" ) )
406  {
407  /* this is the default - use big atom */
408  hmi.lgH2_Chemistry_BigH2 = true;
409  }
410  }
411 
412  /* there is no final branch - if we do not find a keyword, simply
413  * turn on the H2 molecule */
414  return;
415 }
diatomics::lgH2_PAH2_ORNL
bool lgH2_PAH2_ORNL
Definition: h2_priv.h:380
Parser::nMatch
bool nMatch(const char *chKey) const
Definition: parser.h:135
diatomics::n_trace_final
int n_trace_final
Definition: h2_priv.h:402
h2
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
Parser::FFmtRead
double FFmtRead(void)
Definition: parser.cpp:353
diatomics::n_trace_iterations
int n_trace_iterations
Definition: h2_priv.h:403
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
diatomics::lgEnabled
bool lgEnabled
Definition: h2_priv.h:345
diatomics::lgH2_NOISE
bool lgH2_NOISE
Definition: h2_priv.h:383
thirdparty.h
diatomics::coll_source
t_coll_source coll_source[N_X_COLLIDER]
Definition: h2_priv.h:316
diatomics::lgREAD_DATA
bool lgREAD_DATA
Definition: h2_priv.h:252
ParseAtomH2
void ParseAtomH2(Parser &p)
Definition: parse_atom_h2.cpp:13
t_hmi::lgH2_Thermal_BigH2
bool lgH2_Thermal_BigH2
Definition: hmi.h:160
diatomics::xMeanNoise
double xMeanNoise
Definition: h2_priv.h:391
diatomics::n_trace_matrix
int n_trace_matrix
Definition: h2_priv.h:405
diatomics::n_trace_full
int n_trace_full
Definition: h2_priv.h:404
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
Parser::NoNumb
NORETURN void NoNumb(const char *chDesc) const
Definition: parser.cpp:233
Parser
Definition: parser.h:31
diatomics
Definition: h2_priv.h:65
cddefines.h
diatomics::n_elec_states
long int n_elec_states
Definition: h2_priv.h:409
diatomics::xSTDNoise
double xSTDNoise
Definition: h2_priv.h:391
diatomics::nTRACE
int nTRACE
Definition: h2_priv.h:399
diatomics::lgColl_deexec_Calc
bool lgColl_deexec_Calc
Definition: h2_priv.h:359
hmi.h
t_hmi::lgH2_Chemistry_BigH2
bool lgH2_Chemistry_BigH2
Definition: hmi.h:164
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
Parser::lgEOL
bool lgEOL(void) const
Definition: parser.h:98
N_ELEC
const int N_ELEC
Definition: h2_priv.h:27
diatomics::lgH2_grain_deexcitation
bool lgH2_grain_deexcitation
Definition: h2_priv.h:366
diatomics::lgH2_ORH2_ORNL
bool lgH2_ORH2_ORNL
Definition: h2_priv.h:379
diatomics::lgH2_He_ORNL
bool lgH2_He_ORNL
Definition: h2_priv.h:376
h2_priv.h
diatomics::H2_to_H_limit
double H2_to_H_limit
Definition: h2_priv.h:394
parser.h
hmi
t_hmi hmi
Definition: hmi.cpp:5
diatomics::lgH2_ortho_para_coll_on
bool lgH2_ortho_para_coll_on
Definition: h2_priv.h:372
t_coll_source::filename
string filename
Definition: h2_priv.h:62
fixit
void fixit(void)
Definition: service.cpp:991
diatomics::lgColl_dissoc_coll
bool lgColl_dissoc_coll
Definition: h2_priv.h:362
diatomics::lgH2_NOISECOSMIC
bool lgH2_NOISECOSMIC
Definition: h2_priv.h:385
init_genrand
void init_genrand(unsigned long s)
Definition: thirdparty.cpp:2834
diatomics::lgColl_gbar
bool lgColl_gbar
Definition: h2_priv.h:356
h2.h
diatomics::lgLTE
bool lgLTE
Definition: h2_priv.h:369
diatomics::set_numLevelsMatrix
void set_numLevelsMatrix(long numLevels)
Definition: mole_h2_io.cpp:1994
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
deuterium.h
Generated by   doxygen 1.8.17