cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_atom_iso.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2  * others. For conditions of distribution and use see copyright notice in license.txt */
3 /*ParseAtomISO parse information from the atom XX-like command line */
4 #include "cddefines.h"
5 #include "elementnames.h"
6 #include "helike_recom.h"
7 #include "hydrogenic.h"
8 #include "iso.h"
9 #include "parse.h"
10 #include "phycon.h"
11 #include "rfield.h"
12 #include "taulines.h"
13 #include "thirdparty.h"
14 
15 /*ParseAtomISO parse parameters off the XX-like command */
16 void ParseAtomISO(long ipISO, char *chCard )
17 {
18  bool lgEOL;
19  long int i,
20  numLevels,
21  nelem;
22 
23  DEBUG_ENTRY( "ParseAtomISO()" );
24 
25  /* look for the name of an element - if we don't find one do the entire
26  * iso sequence - returns negative number if element not found */
27  nelem = GetElem( chCard );
28 
29  /* H-like Helium is not possible */
30  if( ipISO==ipHE_LIKE && nelem==ipHYDROGEN )
31  {
32  fprintf(ioQQQ," Sorry, H-like He does not exist.\n");
33  cdEXIT(EXIT_FAILURE);
34  }
35 
36  /* collisions - don't pick up the levels collapsed command */
37  if( nMatch("COLL",chCard) && !nMatch("LEVE" , chCard ) )
38  {
39  /* option to turn collisions off, all are
40  * set to 1 in zero. command can accept only one option at a time */
41  if( nMatch("EXCI",chCard) )
42  {
43  /* turn off collisional excitation */
44  iso.lgColl_excite[ipISO] = false;
45  phycon.lgPhysOK = false;
46  }
47  else if( nMatch("IONI",chCard) )
48  {
49  /* turn off collisional ionization */
50  iso.lgColl_ionize[ipISO] = false;
51  phycon.lgPhysOK = false;
52  }
53 
54  else if( nMatch("2S2P",chCard) || ( nMatch("2P2S",chCard) && ipISO == ipH_LIKE ) )
55  {
56  /* >>chng 02 feb 07, change from 2s2p to l-mixing */
57  /* this is the atom h-like collision l-mixing command */
58  fprintf(ioQQQ,"This command changed to ATOM H-LIKE COLLISIONS L-MIXING\n");
59  fprintf(ioQQQ,"I will parse it for now, but may not in the future.\n");
60  /* turn off 2s - 2p collisions */
61  iso.lgColl_l_mixing[ipISO] = false;
62  phycon.lgPhysOK = false;
63  }
64 
65  else if( nMatch("L-MI",chCard) )
66  {
67  if( ipISO == ipH_LIKE )
68  {
69  /* >>chng 02 feb 07, change from 2s2p to l-mixing */
70  /* this is the atom h-like collision l-mixing command */
71  iso.lgColl_l_mixing[ipISO] = false;
72  phycon.lgPhysOK = false;
73  }
74  else if( nMatch("THER",chCard) )
75  {
76  /* use l-mix from
77  *>>refer l-mix all Vrinceanu, D. & Flannery, M. R. 2001, PhysRevA 63, 032701 */
78  if( nMatch("NO T",chCard) )
79  {
80  /* This is the "NO Thermal average" command. It
81  * causes collisions strengths to be evaluated at kT rather than
82  * integrated over a Maxwellian peaked at kT. */
83  iso.lgCS_therm_ave[ipISO] = false;
84  }
85  else
86  {
87  iso.lgCS_therm_ave[ipISO] = true;
88  }
89  }
90  else if( nMatch("PENG",chCard) )
91  {
92  iso.lgCS_Vrinceanu[ipISO] = false;
93  }
94  else if( nMatch(" OFF" , chCard ) )
95  {
96  /* this is the atom xx-like collision l-mixing command */
97  /* turn off same-n collisions */
98  iso.lgColl_l_mixing[ipISO] = false;
99  phycon.lgPhysOK = false;
100  iso.lgCS_Vrinceanu[ipISO] = false;
101  iso.lgCS_therm_ave[ipISO] = false;
102  }
103  else
104  {
105  fprintf( ioQQQ, " needs parameter\n" );
106  cdEXIT(EXIT_FAILURE);
107  }
108  }
109  else if( nMatch(" OFF" , chCard ) )
110  {
111  /* turn everything off, since no keyword given */
112  iso.lgColl_excite[ipISO] = false;
113  iso.lgColl_ionize[ipISO] = false;
114  iso.lgColl_l_mixing[ipISO] = false;
115  phycon.lgPhysOK = false;
116  }
117  else
118  {
119  fprintf( ioQQQ, " needs parameter\n" );
120  cdEXIT(EXIT_FAILURE);
121  }
122  }
123 
124  else if( nMatch("DAMP",chCard) )
125  {
126  if( ipISO == ipHE_LIKE )
127  {
128  fprintf(ioQQQ," Sorry, the DAMPING option is not implemented for the he-like sequence.\n");
129  cdEXIT(EXIT_FAILURE);
130  }
131 
132  /* turn off absorption due to Ly alpha damping wings */
133  hydro.DampOnFac = 0.;
134  }
135 
136  else if( nMatch("DIEL",chCard) )
137  {
138  if( ipISO == ipH_LIKE )
139  {
140  fprintf(ioQQQ," Sorry, but dielectronic recombination onto the h-like sequence is not possible.\n");
141  cdEXIT(EXIT_FAILURE);
142  }
143 
144  /* This sets which set of data to use for dielectronic recombination. */
145  if( nMatch(" OFF",chCard) )
146  {
147  iso.lgDielRecom[ipISO] = false;
148  }
149  else
150  iso.lgDielRecom[ipISO] = true;
151  }
152 
153  else if( nMatch("LEVE",chCard) )
154  {
155  /* the number of levels read in is n, the principal quantum number
156  * only lines with upper levels less than n will be printed */
157 
158  /* number of levels for iso-sequence */
159  /* there are two options here,
160  * when keyword ELEMENT appears, scan off element name and change levels only for
161  * that one.
162  * when there is no ELEMENT then set all in iso to same number */
163 
164  /* lgHydroMalloc is false at start of calculation, set true when space
165  * allocated for the hydrogen and helium lines. Once done we must ignore all
166  * future changes in the number of levels */
167  if( nMatch("PRIN",chCard) )
168  {
169  /* only print - do not change levels */
170  iso.lgPrintNumberOfLevels = true;
171  }
172  else if( !lgHydroMalloc )
173  {
174  i = 5;
175  numLevels = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
176 
177  if( lgEOL )
178  {
179  /* must be a number, or a key, either large (or limit) or small
180  * if no number but LARGER or SMALL then set to preset number */
181  if( nMatch("LARG",chCard) )
182  {
183  /* there is no real limit, but this includes all levels with tabulated rec coefficient */
184  numLevels = RREC_MAXN;
185  }
186 
187  /* this is small or compact keyword */
188  else if( nMatch("SMAL",chCard) || nMatch("COMP",chCard) )
189  {
190  if( ipISO == ipH_LIKE )
191  numLevels = 5;
192  else if( ipISO == ipHE_LIKE )
193  numLevels = 3;
194  }
195 
196  else
197  /* punch out if no number */
198  NoNumb(chCard);
199  }
200  else if( ( numLevels < 3 ) && !nMatch("COLL",chCard) && ipISO == ipH_LIKE )
201  {
202  /* only enforce this limit if not working on collapsed levels */
203  fprintf( ioQQQ, " cannot have fewer than 3 levels, the requested number was %li\n" ,
204  numLevels );
205  fprintf( ioQQQ, " Sorry.\n" );
206  cdEXIT(EXIT_FAILURE);
207  }
208  else if( ( numLevels < 3 ) && !nMatch("COLL",chCard) && ipISO == ipHE_LIKE)
209  {
210  /* only enforce this limit if not working on collapsed levels */
211  fprintf( ioQQQ, " cannot have fewer than 3 levels, the requested number was %li\n" ,
212  numLevels );
213  fprintf( ioQQQ, " Sorry.\n" );
214  cdEXIT(EXIT_FAILURE);
215  }
216 
217  if( ipISO == ipH_LIKE && numLevels > NHYDRO_MAX_LEVEL-2 )
218  {
219  fprintf( ioQQQ, " Not possible to set nhlvl to >NHYDRO_MAX_LEVEL-2= %i\n",
220  NHYDRO_MAX_LEVEL-2 );
221  fprintf( ioQQQ, " change NHYDRO_MAX_LEVEL\n");
222  cdEXIT(EXIT_FAILURE);
223  }
224 
225  /* this is check that alpha transition of highest level is within energy
226  * bounds of continuum */
227  if( ipISO == ipH_LIKE &&
228  ( 2. / POW3((double)numLevels) < rfield.emm ) )
229  {
230  fprintf( ioQQQ, " Not possible to set iso.numLevels_max[ipH_LIKE][ipHYDROGEN] to such a high value, since "
231  "alpha transition not within energy bounds of code\n");
232 
233  fprintf( ioQQQ, " lowest energy is %e and corresponding highest level is %li\n" ,
234  rfield.emm, (long)pow(2./rfield.emm, 0.3333) );
235  cdEXIT(EXIT_FAILURE);
236  }
237 
238  if( nMatch("COLL",chCard) )
239  {
240  if( numLevels < 1 )
241  {
242  fprintf( ioQQQ, "There must be at least one collapsed level.\n");
243  cdEXIT(EXIT_FAILURE);
244  }
245 
246  /* set collapsed levels for entire sequence or just one element */
247  if( nMatch("ELEM",chCard) )
248  {
249  /* check which element is on the line, nelem is atomic number on C scale */
250  nelem = GetElem(chCard);
251  iso.nCollapsed_max[ipISO][nelem] = numLevels;
252  iso_update_num_levels( ipISO, nelem );
253  }
254  else
255  {
256  /* keyword element did not appear, so set all iso species to this number */
257  for( nelem=ipISO; nelem<LIMELM; ++nelem )
258  {
259  iso.nCollapsed_max[ipISO][nelem] = numLevels;
260  iso_update_num_levels( ipISO, nelem );
261  }
262  }
263  }
264 
265  /* we now have the desired number of levels, and it has been fully qualified.
266  * now check whether ELEMENT keyword is on line, if not then set all to
267  * the number, if so then only element in use */
268  else if( nMatch("ELEM",chCard) )
269  {
270  /* check which element is on the line, nelem is atomic number on C scale */
271  nelem = GetElem(chCard);
272  iso.n_HighestResolved_max[ipISO][nelem] = numLevels;
273  iso_update_num_levels( ipISO, nelem );
274  }
275  else
276  {
277  /* keyword element did not appear, so set all species to this number */
278  for( nelem=ipISO; nelem<LIMELM; ++nelem )
279  {
280  iso.n_HighestResolved_max[ipISO][nelem] = numLevels;
281  iso_update_num_levels( ipISO, nelem );
282  }
283  }
284  }
285  }
286 
287  else if( nMatch("ERRO",chCard) && nMatch("GENE" , chCard ) )
288  {
289  /* Rates will be modified by a randomly generated error that falls within
290  * the range specifically set for each rate (or set of rates). */
291  iso.lgRandErrGen[ipISO] = true;
292  i = 5;
293  iso.modelRank[ipISO] = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
294 
295  iso.modelRank[ipISO] = MAX2( 0, iso.modelRank[ipISO] );
296  if( lgEOL )
297  /* Changed to avoid lint complaint */
298  /* iso.modelRank[ipISO] = (unsigned)time(NULL); */
299  iso.modelRank[ipISO] = abs((int)time(NULL));
300 
301  /* this allows a seed that is dependent upon the processor rank
302  * in a parallel run. */
303  /* We add 2 so that the seed is always greater than 1, which would reset the generator. */
304  init_genrand( (unsigned)iso.modelRank[ipISO] + 2);
305  }
306 
307  else if( nMatch(" FSM",chCard) )
308  {
309  if( ipISO == ipH_LIKE )
310  {
311  fprintf(ioQQQ," Sorry, but fine-structure mixing can only be implemented for the He-like sequence.\n");
312  cdEXIT(EXIT_FAILURE);
313  }
314 
315  /* turn on fine structure mixing of spontaneous decays.
316  * >>refer Helike FSM Bauman, R., MacAdams, K., and Ferland, G. (2003). */
317  if( nMatch(" OFF",chCard) )
318  iso.lgFSM[ipISO] = false;
319  else
320  iso.lgFSM[ipISO] = true;
321  }
322 
323  else if( nMatch("GBAR",chCard) )
324  {
325  if( ipISO == ipH_LIKE )
326  {
327  fprintf(ioQQQ," Sorry, the GBAR option is only implemented for the He-like sequence.\n");
328  cdEXIT(EXIT_FAILURE);
329  }
330 
331  /* the HEGBAR command - to change cs of higher levels */
332  /* first turn all off, one will be turned back on */
333  iso.lgCS_Vriens[ipISO] = false;
334  iso.lgCS_None[ipISO] = false;
335  iso.nCS_new[ipISO] = false;
336 
337  /* now turn one on */
338  if( nMatch("VRIE",chCard) )
339  {
340  /* option to change how collisions for unknown levels affect things */
341  iso.lgCS_Vriens[ipISO] = true;
342  }
343  else if( nMatch(" NEW",chCard) )
344  {
345  /* option to change how collisions for unknown levels affect things */
346  i = 5;
347  iso.nCS_new[ipISO] = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
348 
349  /* there are two options, 1 and 2, for which fit - default (no number)
350  * will be 1, the broken power law fit */
351  if( lgEOL )
352  iso.nCS_new[ipISO] = 1;
353 
354  ASSERT( iso.nCS_new[ipISO] );
355  }
356  else if( nMatch(" OFF",chCard) )
357  {
358  /* option to change how collisions for unknown levels affect things */
359  iso.lgCS_None[ipISO] = true;
360  }
361  else
362  {
363  fprintf( ioQQQ, " needs parameter\n" );
364  cdEXIT(EXIT_FAILURE);
365  }
366  }
367 
368 
369  else if( nMatch("LYMA",chCard) )
370  {
371  if( ipISO == ipH_LIKE && nMatch("PUMP",chCard) )
372  {
373  /* >>chng 05 jul 08, separate out Lyman pump commands */
374  if( nMatch(" OFF",chCard) )
375  {
376  /* option to turn off all continuum pumping of Lyman lines */
377  hydro.lgLymanPumping = false;
378  }
379  else if( nMatch("SCALE",chCard) )
380  {
381  /* multiplicative factor for all continuum pumping of H I Lyman lines,
382  * account for possible emission in the line - only affects H I
383  * not entire H-like iso sequence */
384  i = 5;
386  (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
387  /* scale factor is log if negative */
391  }
392  else
393  {
394  fprintf(ioQQQ," Sorry, I didn\'t recognize an option on this ATOM H-LIKE LYMAN PUMP command.\n");
395  fprintf(ioQQQ," The options are \" OFF\", and \"SCALE\".\n");
396  cdEXIT(EXIT_FAILURE);
397  }
398  }
399  else
400  {
401  /* option to set number of "extra" Lyman lines, used for optical depths only */
402  i = 5;
403  iso.nLyman[ipISO] = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
404  if( lgEOL )
405  NoNumb(chCard);
406  }
407  }
408 
409  /* don't interpolate on look-up tables but compute recombination on the fly instead */
410  else if( nMatch("RECO",chCard) &&
411  nMatch(" NO ",chCard) && nMatch("INTE",chCard) )
412  {
413  /* flag set by atom xx-like no recombination interp command,
414  * says to generate recombination coefficients
415  * on the fly */
416  iso.lgNoRecombInterp[ipISO] = true;
417  }
418 
419  else if( nMatch("REDI",chCard) )
420  {
421  int ipRedis=0;
422  /* there are three functions, PRD_, CRD_, and CRDW,
423  * representing partial redistribution,
424  * complete redistribution (doppler core only, no wings)
425  * and complete with wings */
426  /* partial redistribution */
427  if( nMatch(" PRD",chCard) )
428  {
429  ipRedis = ipPRD;
430  }
431  /* complete redistribution no wings */
432  else if( nMatch(" CRD",chCard) )
433  {
434  ipRedis = ipCRD;
435  }
436  /* complete redistribution with wings */
437  else if( nMatch("CRDW",chCard) )
438  {
439  ipRedis = ipCRDW;
440  }
441 
442  /* if not SHOW option (handled below) then we have a problem */
443  else if( !nMatch("SHOW",chCard) )
444  {
445  fprintf(ioQQQ," There should have been a second keyword on this command.\n");
446  fprintf(ioQQQ," Options are _PRD, _CRD, CRDW (_ is space). Sorry.\n");
447  cdEXIT(EXIT_FAILURE);
448  }
449 
450  /* resonance lines - not Lya*/
451  if( nMatch("ALPH",chCard) )
452  {
453  iso.ipLyaRedist[ipISO] = ipRedis;
454  }
455  /* Lya itself */
456  if( nMatch("RESO",chCard) )
457  {
458  iso.ipResoRedist[ipISO] = ipRedis;
459  }
460  /* subordinate lines */
461  else if( nMatch("SUBO",chCard) )
462  {
463  iso.ipSubRedist[ipISO] = ipRedis;
464  }
465  /* the show option, say what we are assuming */
466  else if( nMatch("SHOW",chCard) )
467  {
468  fprintf(ioQQQ," Ly a is ");
469  if( iso.ipLyaRedist[ipISO] ==ipCRDW )
470  {
471  fprintf(ioQQQ,"complete redistribution with wings\n");
472  }
473  else if( iso.ipLyaRedist[ipISO] ==ipCRD )
474  {
475  fprintf(ioQQQ,"complete redistribution with core only.\n");
476  }
477  else if( iso.ipLyaRedist[ipISO] ==ipPRD )
478  {
479  fprintf(ioQQQ,"partial redistribution.\n");
480  }
481  else if( iso.ipLyaRedist[ipISO] ==ipLY_A )
482  {
483  fprintf(ioQQQ,"special Lya.\n");
484  }
485  else
486  {
487  fprintf(ioQQQ," PROBLEM Impossible value for iso.ipLyaRedist.\n");
488  TotalInsanity();
489  }
490 
491  fprintf(ioQQQ," Other %s resonance lines are ",
492  elementnames.chElementSym[ipISO] );
493 
494  if( iso.ipResoRedist[ipISO] ==ipCRDW )
495  {
496  fprintf(ioQQQ,"complete redistribution with wings\n");
497  }
498  else if( iso.ipResoRedist[ipISO] ==ipCRD )
499  {
500  fprintf(ioQQQ,"complete redistribution with core only.\n");
501  }
502  else if( iso.ipResoRedist[ipISO] ==ipPRD )
503  {
504  fprintf(ioQQQ,"partial redistribution.\n");
505  }
506  else
507  {
508  fprintf(ioQQQ," PROBLEM Impossible value for iso.ipResoRedist.\n");
509  TotalInsanity();
510  }
511 
512  fprintf(ioQQQ," %s subordinate lines are ",
513  elementnames.chElementSym[ipISO] );
514 
515  if( iso.ipSubRedist[ipISO] ==ipCRDW )
516  {
517  fprintf(ioQQQ,"complete redistribution with wings\n");
518  }
519  else if( iso.ipSubRedist[ipISO] ==ipCRD )
520  {
521  fprintf(ioQQQ,"complete redistribution with core only.\n");
522  }
523  else if( iso.ipSubRedist[ipISO] ==ipPRD )
524  {
525  fprintf(ioQQQ,"partial redistribution.\n");
526  }
527  else
528  {
529  fprintf(ioQQQ," PROBLEM Impossible value for iso.ipSubRedist.\n");
530  TotalInsanity();
531  }
532  }
533  else
534  {
535  fprintf(ioQQQ," here should have been another keyword on this command.\n");
536  fprintf(ioQQQ," Options are ALPHA, RESONANCE, SUBORDINATE. Sorry.\n");
537  cdEXIT(EXIT_FAILURE);
538  }
539  }
540 
541  /* which type of matrix solution, populations or lowt */
542  else if( nMatch("MATR",chCard) )
543  {
544  if( nMatch(" POP",chCard) )
545  {
546  strcpy( iso.chTypeAtomSet[ipISO] , "POPU" );
547  }
548  else if( nMatch(" LOW",chCard) )
549  {
550  strcpy( iso.chTypeAtomSet[ipISO] , "LOWT" );
551  }
552  else
553  {
554  fprintf( ioQQQ, " There should have been a keyword on the MATRIX option.\n");
555  fprintf( ioQQQ, " The keywords are POPulations, and LOWte.\n");
556  cdEXIT(EXIT_FAILURE);
557  }
558  }
559 
560  else if( nMatch("TOPO",chCard) )
561  {
562  if( ipISO == ipH_LIKE )
563  {
564  fprintf( ioQQQ, " The behavior of the TOPOFF option for the H-like sequence has changed since version 07.02.\n");
565  fprintf( ioQQQ, " Originally, it accepted a number and distributed \"topoff\" equally into that many levels at the top of the model.\n");
566  fprintf( ioQQQ, " Now, the option works the same as in the He-like sequence, which is to turn off topoff altogether.\n" );
567  }
568 
569  if( nMatch(" OFF",chCard) )
570  {
571  iso.lgTopoff[ipISO] = false;
572  fprintf( ioQQQ, "ISO %li TOPOFF is OFF\n", ipISO );
573  }
574  else
575  iso.lgTopoff[ipISO] = true;
576  }
577 
578 
579  else
580  {
581  fprintf( ioQQQ, " There should have been a keyword - STOP\n" );
582  cdEXIT(EXIT_FAILURE);
583  }
584  return;
585 }

Generated for cloudy by doxygen 1.8.4