cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
parse_element.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 /*ParseElement parse options on element command */
4 #include "cddefines.h"
5 #include "optimize.h"
6 #include "abund.h"
7 #include "dense.h"
8 #include "input.h"
9 #include "iso.h"
10 #include "hmi.h"
11 #include "mole.h"
12 #include "elementnames.h"
13 #include "parse.h"
14 
15 void ParseElement(char *chCard )
16 {
17  /* this will be used to remember how many levels are active in any element we turn off,
18  * so that we can retain state if turned back on */
19  static bool lgFirst = true;
20  static long levels[NISO][LIMELM];
21  char chCap[INPUT_LINE_LENGTH];
22  bool lgEOL,
23  lgEnd,
24  lgHIT;
25 
26  bool lgForceLog=false, lgForceLinear=false;
27 
28  long int i,
29  nelem,
30  j,
31  k;
32  double param;
33 
34  DEBUG_ENTRY( "ParseElement()" );
35 
36  /* zero out this array if first call */
37  if( lgFirst )
38  {
39  lgFirst = false;
40  for( i=0; i<NISO; ++i )
41  {
42  for( nelem=0; nelem<LIMELM; ++nelem )
43  {
44  levels[i][nelem] = 0;
45  }
46  }
47  }
48  /* say that abundances have been changed */
49  abund.lgAbnSolar = false;
50 
51  /* read in list of elements for abundances command */
52  if( nMatch("READ",chCard) )
53  {
54  abund.npSolar = 0;
55  for( i=1; i < LIMELM; i++ )
56  {
57  input_readarray(chCard,&lgEnd);
58 
59  if( lgEnd )
60  {
61  fprintf( ioQQQ, " Hit EOF while reading element list; use END to end list.\n" );
62  cdEXIT(EXIT_FAILURE);
63  }
64 
65  strcpy( chCap, chCard );
66  caps(chCap);
67 
68  /* this would be a line starting with END to say end on list */
69  if( strncmp(chCap,"END",3) == 0 )
70  {
71  return;
72  }
73 
74  j = 1;
75  lgHIT = false;
76  while( j < LIMELM && !lgHIT )
77  {
78  j += 1;
79 
80  if( strncmp( chCap,elementnames.chElementNameShort[j-1] , 4) == 0 )
81  {
82  abund.npSolar += 1;
83  abund.ipSolar[abund.npSolar-1] = j;
84  lgHIT = true;
85  }
86  }
87 
88  if( !lgHIT )
89  {
90  fprintf( ioQQQ, "%80.80s\n", chCard );
91  fprintf( ioQQQ, " Sorry, but I did not recognize element name on this line.\n" );
92  fprintf( ioQQQ, " Here is the list of names I recognize.\n" );
93  fprintf( ioQQQ, " " );
94 
95  for( k=2; k <= LIMELM; k++ )
96  {
97  fprintf( ioQQQ, "%4.4s\n", elementnames.chElementNameShort[k-1] );
98  }
99 
100  cdEXIT(EXIT_FAILURE);
101  }
102  }
103 
104  /* fell through, make sure one more line, the end line */
105  input_readarray(chCard,&lgEnd);
106  strcpy( chCap, chCard );
107  caps(chCap);
108 
109  if( strncmp(chCap,"END",3) == 0 )
110  {
111  return;
112  }
113 
114  else
115  {
116  fprintf( ioQQQ, " Too many elements were entered.\n" );
117  fprintf( ioQQQ, " I only know about%3d elements.\n",
118  LIMELM );
119  fprintf( ioQQQ, " Sorry.\n" );
120  cdEXIT(EXIT_FAILURE);
121  }
122  }
123 
124  /* find which element - will be used in remainder of routine to
125  * adjust aspects of this element */
126  nelem = GetElem(chCard );
127 
128  bool lgElementSet = false;
129  if( nelem >= ipHYDROGEN )
130  {
131  lgElementSet = true;
132  /* nelem is now the atomic number of the element, and must be correct */
133  ASSERT( nelem>=0 && nelem < LIMELM );
134  }
135 
136  /* look for log or linear keywords */
137  if( nMatch(" LOG",chCard) )
138  lgForceLog = true;
139  else if( nMatch("LINE",chCard) )
140  lgForceLinear = true;
141 
142  i = 4;
143  param = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
144 
145  if( nMatch("SCAL",chCard) )
146  {
147  /* enter abundance as scale factor, relative to what is in now */
148  if( lgEOL )
149  {
150  fprintf( ioQQQ, " There must be a number on this line.\n" );
151  fprintf( ioQQQ, " Sorry.\n" );
152  cdEXIT(EXIT_FAILURE);
153  }
154 
155  else
156  {
157 
158  if( !lgElementSet )
159  {
160  fprintf( ioQQQ,
161  " ParseElement did not find an element on the following line:\n" );
162  fprintf( ioQQQ,
163  "%80.80s\n", chCard );
164  cdEXIT(EXIT_FAILURE);
165  }
166  /* interpret as log unless forced linear */
167  if( (lgForceLog || param <= 0.) && !lgForceLinear )
168  {
169  /* option for log of scale factor */
170  param = pow(10.,param);
171  }
172  abund.ScaleElement[nelem] = (realnum)param;
173  }
174  }
175 
176  else if( nMatch("ABUN",chCard) )
177  {
178  /* log of actual abundance */
179  if( lgEOL )
180  {
181  fprintf( ioQQQ, " There must be a number on this line.\n" );
182  fprintf( ioQQQ, " Sorry.\n" );
183  cdEXIT(EXIT_FAILURE);
184  }
185 
186  else
187  {
188 
189  if( !lgElementSet )
190  {
191  fprintf( ioQQQ,
192  " ParseElement did not find an element on the following line:\n" );
193  fprintf( ioQQQ,
194  "%80.80s\n", chCard );
195  cdEXIT(EXIT_FAILURE);
196  }
197  if( lgForceLinear )
198  {
199  abund.solar[nelem] = (realnum)param;
200  }
201  else
202  {
203  abund.solar[nelem] = (realnum)pow(10.,param);
204  }
205 
206  if( abund.solar[nelem]> 1. )
207  {
208  fprintf( ioQQQ,
209  " Please check the abundance of this element. It seems high to me.\n" );
210  }
211  }
212  }
213 
214  else if( nMatch(" OFF",chCard) )
215  {
216  /* keyword LIMIT sets lower limit to abundance of element
217  * that will be included */
218  /* option to turn off this element, set abundance to zero */
219  if( nMatch( "LIMI",chCard) )
220  {
221  if( lgEOL )
222  {
223  fprintf( ioQQQ, " There must be an abundances on the ELEMENT OFF LIMIT command.\n" );
224  fprintf( ioQQQ, " Sorry.\n" );
225  cdEXIT(EXIT_FAILURE);
226  }
227  dense.AbundanceLimit = (realnum)pow(10., param);
228  }
229  else
230  {
231  if( !lgElementSet )
232  {
233  fprintf( ioQQQ,
234  " ParseElement did not find an element on the following line:\n" );
235  fprintf( ioQQQ,
236  "%80.80s\n", chCard );
237  cdEXIT(EXIT_FAILURE);
238  }
239  /* specify element name */
240  dense.lgElmtOn[nelem] = false;
241  /* another flag that element is off */
242  dense.IonHigh[nelem] = -1;
243  dense.lgSetIoniz[nelem] = false;
244  /* set no levels for all elements that are turned off, except for helium itself, which always exists */
245  if( nelem > ipHELIUM )
246  {
247  levels[ipHYDROGEN][nelem] = iso.numLevels_max[ipH_LIKE][nelem];
248  levels[ipHELIUM][nelem] = iso.numLevels_max[ipHE_LIKE][nelem];
249  iso.numLevels_max[ipH_LIKE][nelem] = 0;
250  iso.numLevels_max[ipHE_LIKE][nelem] = 0;
251  }
252 
253  if( nelem == ipHYDROGEN )
254  {
255  fprintf( ioQQQ, " It is not possible to turn hydrogen off.\n" );
256  fprintf( ioQQQ, " Sorry.\n" );
257  cdEXIT(EXIT_FAILURE);
258  }
259  }
260  }
261 
262  /* specify an ionization distribution */
263  else if( nMatch("IONI",chCard) )
264  {
265  bool lgLogSet = false;
266  int ion;
267  int ihi , low;
268 
269 
270  if( !lgElementSet )
271  {
272  fprintf( ioQQQ,
273  " ParseElement did not find an element on the following line:\n" );
274  fprintf( ioQQQ,
275  "%80.80s\n", chCard );
276  cdEXIT(EXIT_FAILURE);
277  }
278  if( !dense.lgElmtOn[nelem] )
279  {
280  fprintf(ioQQQ,"Sorry, you cannot set the ionization of %s since it has been turned off.\n" ,
281  elementnames.chElementName[nelem] );
282  cdEXIT(EXIT_FAILURE);
283  }
284 
285  i = 4;
286  /* option to specify ionization distribution for this element */
287  dense.lgSetIoniz[nelem] = true;
288  ion = 0;
289  while( ion<nelem+2 )
290  {
291  /* the ionization fractions that are set when above set true,
292  * gas phase abundance is this times total abundance
293  * Ionization fraction for [nelem][ion] */
294  dense.SetIoniz[nelem][ion] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
295 
296  if( lgEOL )
297  break;
298 
299  /* all are log if any are negative */
300  if( dense.SetIoniz[nelem][ion] < 0. )
301  lgLogSet = true;
302  ++ion;
303  }
304 
305  /* zero out ones not specified */
306  for( i=ion; i<nelem+2; ++i )
307  dense.SetIoniz[nelem][i] = 0.;
308 
309  /* convert rest to linear if any were negative */
310  if( lgLogSet )
311  {
312  for( i=0; i<ion; ++i )
313  dense.SetIoniz[nelem][i] = (realnum)pow((realnum)10.f, dense.SetIoniz[nelem][i]);
314  }
315 
316  /* now check that no zero abundances exist between lowest and highest non-zero values */
317  low = 0;
318  while( dense.SetIoniz[nelem][low]==0 && low < nelem+1 )
319  ++low;
320  ihi = nelem+1;
321  while( dense.SetIoniz[nelem][ihi]==0 && ihi > low )
322  --ihi;
323 
324  if( ihi==low && dense.SetIoniz[nelem][ihi]==0 )
325  {
326  fprintf(ioQQQ," element ionization command has all zero ionization fractions. This is not possible.\n Sorry\n");
327  cdEXIT(EXIT_FAILURE);
328  }
329  for(i=low; i<=ihi; ++i )
330  {
331  if( dense.SetIoniz[nelem][i]==0 )
332  {
333  fprintf(ioQQQ," element abundance command has zero abundance between positive values. This is not possible.\n Sorry\n");
334  cdEXIT(EXIT_FAILURE);
335  }
336  }
337  /* >>chng 05 mar 13, add this
338  * finally turn off any part of the chemical networks that might use this element,
339  * since chemistry must come into equilibrium with the ionization, and vice vese -
340  * ion distribution cannot come into equilibrium with chemistry when ions are set,
341  * so much turn off chemistry */
342  if( nelem==ipHYDROGEN && !hmi.lgNoH2Mole )
343  {
344  fprintf(ioQQQ," Hydrogen ionization has been set, H chemistry is disabled.\n");
345  /* turn off only H2 */
346  hmi.lgNoH2Mole = true;
347  }
348  /* loop over all atoms in co chem net - these occur as the last NUM_ELEMENTS
349  * elements in the full set of mole.num_heavy_molec */
350  /* >>chng 06 mar 19, do not print if network already turned off */
351  if( mole.lgElem_in_chemistry[nelem] && !co.lgNoCOMole )
352  {
353  fprintf(ioQQQ," The ionization of an element in the CO chemistry network has been set, CO chemistry is disabled.\n");
354  /* turn off CO network */
355  co.lgNoCOMole = true;
356  }
357  }
358 
359  /* option to turn element on - some ini files turn off elements,
360  * this can turn them back on */
361  else if( nMatch(" ON ",chCard) )
362  {
363 
364  if( !lgElementSet )
365  {
366  fprintf( ioQQQ,
367  " ParseElement did not find an element on the following line:\n" );
368  fprintf( ioQQQ,
369  "%80.80s\n", chCard );
370  cdEXIT(EXIT_FAILURE);
371  }
372  /* option to turn off this element, set abundance to zero */
373  dense.lgElmtOn[nelem] = true;
374  /* reset levels to default if they were ever turned off with element off command */
375  if( levels[ipHYDROGEN][nelem] )
376  {
377  iso.numLevels_max[ipH_LIKE][nelem] = levels[ipHYDROGEN][nelem];
378  iso.numLevels_max[ipHE_LIKE][nelem] = levels[ipHELIUM][nelem];
379  }
380  }
381 
382  else if( nMatch("TABL",chCard) )
383  {
384 
385  if( !lgElementSet )
386  {
387  fprintf( ioQQQ,
388  " ParseElement did not find an element on the following line:\n" );
389  fprintf( ioQQQ,
390  "%80.80s\n", chCard );
391  cdEXIT(EXIT_FAILURE);
392  }
393  /* read in table of position-dependent abundances for a particular element. */
394  abund.lgAbunTabl[nelem] = true;
395 
396  /* general flag saying this option turned on */
397  abund.lgAbTaON = true;
398 
399  /* does the table give depth or radius? keyword DEPTH
400  * says it is depth, default is radius */
401  if( nMatch("DEPT",chCard) )
402  abund.lgAbTaDepth[nelem] = true;
403  else
404  abund.lgAbTaDepth[nelem] = false;
405 
406  /* make sure not trying to change hydrogen */
407  if( nelem == ipHYDROGEN )
408  {
409  fprintf( ioQQQ, " cannot change abundance of hydrogen.\n" );
410  fprintf( ioQQQ, " Sorry.\n" );
411  cdEXIT(EXIT_FAILURE);
412  }
413 
414  /* read pair giving radius and abundance */
415  input_readarray(chCard,&lgEnd);
416  i = 1;
417  abund.AbTabRad[0][nelem] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
418  abund.AbTabFac[0][nelem] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
419 
420  if( lgEOL )
421  {
422  fprintf( ioQQQ, " no pairs entered - cannot interpolate\n" );
423  cdEXIT(EXIT_FAILURE);
424  }
425 
426  /* number of points in the table */
427  abund.nAbunTabl = 2;
428 
429  lgEnd = false;
430  /* LIMTAB is limit to number of points we can store and is
431  * set to 500 in abundances */
432  while( !lgEnd && abund.nAbunTabl < LIMTABD )
433  {
434  input_readarray(chCard,&lgEnd);
435  if( !lgEnd )
436  {
437  /* convert first 4 char to caps, into chCap */
438  cap4(chCap , chCard);
439  if( strncmp(chCap,"END",3) == 0 )
440  lgEnd = true;
441  }
442 
443  /* lgEnd may have been set within above if, if end line encountered*/
444  if( !lgEnd )
445  {
446  i = 1;
447  abund.AbTabRad[abund.nAbunTabl-1][nelem] =
448  (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
449 
450  abund.AbTabFac[abund.nAbunTabl-1][nelem] =
451  (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
452  abund.nAbunTabl += 1;
453  }
454  }
455 
456  abund.nAbunTabl -= 1;
457 
458  /* now check that radii are in increasing order */
459  for( i=1; i < abund.nAbunTabl; i++ )
460  {
461  /* the radius values are assumed to be strictly increasing */
462  if( abund.AbTabRad[i][nelem] <= abund.AbTabRad[i-1][nelem] )
463  {
464  fprintf( ioQQQ, "ParseElement: TABLE ELEMENT TABLE radii "
465  "must be in increasing order\n" );
466  cdEXIT(EXIT_FAILURE);
467  }
468  }
469  }
470 
471  else
472  {
473  fprintf( ioQQQ, "ParseElement: ELEMENT command - there must be "
474  "a keyword on this line.\n" );
475  fprintf( ioQQQ, " The keys I know about are TABLE, SCALE, _OFF, "
476  "_ON_, IONIZATION, and ABUNDANCE.\n" );
477  fprintf( ioQQQ, " Sorry.\n" );
478  cdEXIT(EXIT_FAILURE);
479  }
480 
481  /* vary option */
482  if( optimize.lgVarOn )
483  {
485  /* pointer to where to write */
487 
488  if( nMatch("SCAL",chCard) )
489  {
490  /* vary scale factor */
491  sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s SCALE %%f LOG",
493 
494  /* param is linear scale factor */
495  optimize.vparm[0][optimize.nparm] = (realnum)log10(param);
496  optimize.vincr[optimize.nparm] = 0.2f;
497  }
498 
499  else if( nMatch("ABUN",chCard) )
500  {
501  /* vary absolute abundance */
502  sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s ABUND %%f LOG ",
504 
505  /* param is log of abundance by number relative to hydrogen */
506  optimize.vparm[0][optimize.nparm] = (realnum)param;
507  optimize.vincr[optimize.nparm] = 0.2f;
508  }
509  ++optimize.nparm;
510  }
511  return;
512 }

Generated for cloudy by doxygen 1.8.1.2