cloudy  trunk
 All Data Structures Namespaces 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-2022 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 "called.h"
10 #include "iso.h"
11 #include "mole.h"
12 #include "elementnames.h"
13 #include "parser.h"
14 #include "deuterium.h"
15 
16 
18 {
19  /* this will be used to remember how many levels are active in any element we turn off,
20  * so that we can retain state if turned back on */
21  static bool lgFirst = true;
22  static long levels[NISO][LIMELM];
23  bool lgEnd,
24  lgHIT;
25 
26  bool lgForceLog=false, lgForceLinear=false;
27 
28  long int nelem,
29  j,
30  k;
31  double param=0.0;
32 
33  DEBUG_ENTRY( "ParseElement()" );
34 
35  /* zero out this array if first call */
36  if( lgFirst )
37  {
38  lgFirst = false;
39  for( long i=0; i<NISO; ++i )
40  {
41  for( nelem=0; nelem<LIMELM; ++nelem )
42  {
43  levels[i][nelem] = 0;
44  }
45  }
46  }
47  /* say that abundances have been changed */
48  abund.lgAbnSolar = false;
49 
50  /* read in list of elements for abundances command */
51  if( p.nMatch("READ") )
52  {
53  abund.npSolar = 0;
54  for( long i=1; i < LIMELM; i++ )
55  {
56  p.getline();
57 
58  if( p.m_lgEOF )
59  {
60  fprintf( ioQQQ, " Hit EOF while reading element list; use END to end list.\n" );
62  }
63 
64  /* this would be a line starting with END to say end on list */
65  if( p.hasCommand("END") )
66  {
67  return;
68  }
69 
70  j = 1;
71  lgHIT = false;
72  while( j < LIMELM && !lgHIT )
73  {
74  j += 1;
75 
77  {
78  abund.npSolar += 1;
79  abund.ipSolar[abund.npSolar-1] = j;
80  lgHIT = true;
81  }
82  }
83 
84  if( !lgHIT )
85  {
86  p.PrintLine( ioQQQ);
87  fprintf( ioQQQ, " Sorry, but I did not recognize element name on this line.\n" );
88  fprintf( ioQQQ, " Here is the list of names I recognize.\n" );
89  fprintf( ioQQQ, " " );
90 
91  for( k=2; k <= LIMELM; k++ )
92  {
93  fprintf( ioQQQ, "%4.4s\n", elementnames.chElementNameShort[k-1] );
94  }
95 
97  }
98  }
99 
100  /* fell through, make sure one more line, the end line */
101  p.getline();
102 
103  if( p.hasCommand("END") )
104  {
105  return;
106  }
107 
108  else
109  {
110  fprintf( ioQQQ, " Too many elements were entered.\n" );
111  fprintf( ioQQQ, " I only know about%3d elements.\n",
112  LIMELM );
113  fprintf( ioQQQ, " Sorry.\n" );
115  }
116  }
117 
118  /* find which element - will be used in remainder of routine to
119  * adjust aspects of this element */
120  nelem = p.GetElem( );
121 
122  bool lgElementSet = false;
123  if( nelem >= ipHYDROGEN )
124  {
125  lgElementSet = true;
126  /* nelem is now the atomic number of the element, and must be correct */
127  ASSERT( nelem>=0 && nelem < LIMELM );
128  }
129 
130  /* look for log or linear keywords */
131  if( p.nMatch(" LOG") )
132  lgForceLog = true;
133  else if( p.nMatch("LINE") )
134  lgForceLinear = true;
135 
136  if( p.nMatch("SCAL") )
137  {
138  param = p.FFmtRead();
139 
140  /* enter abundance as scale factor, relative to what is in now */
141  if( p.lgEOL() )
142  {
143  fprintf( ioQQQ, " There must be a number on this line.\n" );
144  fprintf( ioQQQ, " Sorry.\n" );
146  }
147 
148  else
149  {
150  if( !lgElementSet )
151  {
152  fprintf( ioQQQ,
153  " ParseElement did not find an element on the following line:\n" );
154  p.PrintLine( ioQQQ );
156  }
157  /* interpret as log unless forced linear */
158  if( (lgForceLog || param <= 0.) && !lgForceLinear )
159  {
160  /* option for log of scale factor */
161  param = exp10(param);
162  }
163  abund.ScaleElement[nelem] = (realnum)param;
164  }
165  }
166 
167  else if( p.nMatch("ABUN") )
168  {
169  param = p.FFmtRead();
170 
171  /* log of actual abundance */
172  if( p.lgEOL() )
173  {
174  fprintf( ioQQQ, " There must be a number on this line.\n" );
175  fprintf( ioQQQ, " Sorry.\n" );
177  }
178 
179  else
180  {
181  if( !lgElementSet )
182  {
183  fprintf( ioQQQ,
184  " ParseElement did not find an element on the following line:\n" );
185  p.PrintLine( ioQQQ );
187  }
188  /* interpret as log unless forced linear */
189  if( !lgForceLinear )
190  {
191  /* option for log of scale factor */
192  if( param > log10(BIGFLOAT) )
193  {
194  fprintf(ioQQQ,"log of abundance was entered, but it too large for this cpu.\n");
196  }
197  param = exp10(param);
198  }
199  abund.solar[nelem] = (realnum)param;
200 
201  if( abund.solar[nelem] > 1. && called.lgTalk )
202  {
203  fprintf( ioQQQ,
204  " Please check the abundance of this element. It seems high to me.\n" );
205  }
206  }
207  }
208 
209  else if( p.nMatch("ISOT") )
210  {
211  if( p.nMatch(" ALL") )
212  {
213  for( int i = ipHELIUM; i < LIMELM; i++ )
214  mole_global.lgTreatIsotopes[i] = true;
215  }
216  else
217  {
218  if( !lgElementSet )
219  {
220  fprintf( ioQQQ,
221  " ParseElement did not find an element on the following line:\n" );
222  p.PrintLine( ioQQQ );
224  }
225 
226  if( !dense.lgElmtOn[nelem] )
227  {
228  fprintf(ioQQQ,"Sorry, you cannot set the isotope fractions of %s since it has been turned off.\n" ,
229  elementnames.chElementName[nelem] );
231  }
232 
233  /* option to specify ionization distribution for this element */
234  int iso = 0;
235  while( true )
236  {
237  int Aiso = (int)p.FFmtRead();
238  if( p.lgEOL() ) break;
239 
240  realnum Fiso = (realnum)p.FFmtRead();
241  if( p.lgEOL() )
242  p.NoNumb( "isotope abundance" );
243 
244 
245  int j = abund.IsoAbn[nelem].setAbn( Aiso, Fiso );
246  if( j == -1 )
247  {
248  fprintf(ioQQQ, "Element Isotope: For %s, illegal isotope mass number: %d\n",
249  elementnames.chElementName[nelem], Aiso);
250  fprintf(ioQQQ, "Choose one of:\t");
251  abund.IsoAbn[nelem].prtIso( ioQQQ );
253  }
254 
255  iso++;
256  }
257 
258  if( iso != abund.IsoAbn[nelem].getNiso() )
259  {
260  fprintf(ioQQQ, "Element Isotope: %s requires %d isotope pairs to be specified, but found %d\n",
261  elementnames.chElementName[nelem], abund.IsoAbn[nelem].getNiso(), iso);
262  fprintf(ioQQQ, "Enter fractionations for all isotopes, eg:\n");
263  fprintf(ioQQQ, "element %s isotopes", elementnames.chElementName[nelem]);
264  abund.IsoAbn[nelem].prtIsoPairs( ioQQQ );
266  }
267 
268  if( abund.IsoAbn[nelem].isAnyIllegal() )
269  {
270  fprintf(ioQQQ, "Element Isotope: Non-positive isotope fractions are illegal.\n");
271  fprintf(ioQQQ, "Read: %s\t iso:", elementnames.chElementName[nelem]);
272  abund.IsoAbn[nelem].prtIsoPairs( ioQQQ );
274  }
275 
276  abund.IsoAbn[nelem].normAbn( );
277  // abund.IsoAbn[nelem].prtIsoPairs( stdout );
278 
279  mole_global.lgTreatIsotopes[nelem] = true;
280  if( nelem == ipHYDROGEN )
281  {
282  deut.lgElmtOn = true;
283  if( abund.IsoAbn[nelem].getAbn(2) == 0.0 )
284  deut.lgElmtOn = false;
285  }
286  }
287  }
288 
289  else if( p.nMatch(" OFF") )
290  {
291  /* keyword LIMIT sets lower limit to abundance of element
292  * that will be included */
293  /* option to turn off this element, set abundance to zero */
294  if( p.nMatch( "LIMI") )
295  {
296  param = p.FFmtRead();
297  if( p.lgEOL() )
298  {
299  fprintf( ioQQQ, " There must be an abundances on the ELEMENT OFF LIMIT command.\n" );
300  fprintf( ioQQQ, " Sorry.\n" );
302  }
303  dense.AbundanceLimit = (realnum)exp10(param);
304  }
305  else
306  {
307  if( !lgElementSet )
308  {
309  fprintf( ioQQQ,
310  " ParseElement did not find an element on the following line:\n" );
311  p.PrintLine ( ioQQQ );
313  }
314  /* flags saying element off, explicitly set by user */
315  dense.lgElmtOn[nelem] = false;
316  dense.lgElmtSetOff[nelem] = true;
317  /* another flag that element is off */
318  dense.IonHigh[nelem] = -1;
319  dense.lgSetIoniz[nelem] = false;
320  /* set no levels for all elements that are turned off, except for helium itself, which always exists */
321  if( nelem > ipHELIUM )
322  {
323  levels[ipHYDROGEN][nelem] = iso_sp[ipH_LIKE][nelem].numLevels_max;
324  levels[ipHELIUM][nelem] = iso_sp[ipHE_LIKE][nelem].numLevels_max;
325  iso_sp[ipH_LIKE][nelem].numLevels_max = 0;
326  iso_sp[ipHE_LIKE][nelem].numLevels_max = 0;
327  }
328 
329  if( nelem == ipHYDROGEN )
330  {
331  fprintf( ioQQQ, " It is not possible to turn hydrogen off.\n" );
332  fprintf( ioQQQ, " Sorry.\n" );
334  }
335  }
336  }
337 
338  /* specify an ionization distribution */
339  else if( p.nMatch("IONI") )
340  {
341  bool lgLogSet = false;
342  int ion;
343  int ihi , low;
344 
345 
346  if( !lgElementSet )
347  {
348  fprintf( ioQQQ,
349  " ParseElement did not find an element on the following line:\n" );
350  p.PrintLine( ioQQQ );
352  }
353  if( !dense.lgElmtOn[nelem] )
354  {
355  fprintf(ioQQQ,"Sorry, you cannot set the ionization of %s since it has been turned off.\n" ,
356  elementnames.chElementName[nelem] );
358  }
359 
360  /* option to specify ionization distribution for this element */
361  dense.lgSetIoniz[nelem] = true;
362  ion = 0;
363  while( ion<nelem+2 )
364  {
365  /* the ionization fractions that are set when above set true,
366  * gas phase abundance is this times total abundance
367  * Ionization fraction for [nelem][ion] */
368  dense.SetIoniz[nelem][ion] = (realnum)p.FFmtRead();
369 
370  if( p.lgEOL() )
371  break;
372 
373  /* all are log if any are negative */
374  if( dense.SetIoniz[nelem][ion] < 0. )
375  lgLogSet = true;
376  ++ion;
377  }
378  param = dense.SetIoniz[nelem][0];
379 
380  /* zero out ones not specified */
381  for( long i=ion; i<nelem+2; ++i )
382  dense.SetIoniz[nelem][i] = 0.;
383 
384  /* convert rest to linear if any were negative */
385  if( lgLogSet )
386  {
387  for( long i=0; i<ion; ++i )
388  dense.SetIoniz[nelem][i] = exp10(dense.SetIoniz[nelem][i]);
389  }
390 
391  /* now check that no zero abundances exist between lowest and highest non-zero values */
392  low = 0;
393  while( dense.SetIoniz[nelem][low]==0 && low < nelem+1 )
394  ++low;
395  ihi = nelem+1;
396  while( dense.SetIoniz[nelem][ihi]==0 && ihi > low )
397  --ihi;
398 
399  if( ihi==low && dense.SetIoniz[nelem][ihi]==0 )
400  {
401  fprintf(ioQQQ," element ionization command has all zero ionization fractions. This is not possible.\n Sorry\n");
403  }
404  realnum AbundSum=0.;
405  for(long i=low; i<=ihi; ++i )
406  {
407  if( dense.SetIoniz[nelem][i]==0 )
408  {
409  fprintf(ioQQQ," element abundance command has zero abundance between positive values. This is not possible.\n Sorry\n");
411  }
412  AbundSum += dense.SetIoniz[nelem][i];
413  }
414 
415  // renormalize so that sum is unity - needed for various sum rules
416  for(long i=low; i<=ihi; ++i )
417  dense.SetIoniz[nelem][i] /= AbundSum;
418 
419  }
420 
421  /* option to turn element on - some ini files turn off elements,
422  * this can turn them back on */
423  else if( p.nMatch(" ON ") )
424  {
425 
426  if( !lgElementSet )
427  {
428  fprintf( ioQQQ,
429  " ParseElement did not find an element on the following line:\n" );
430  p.PrintLine( ioQQQ );
432  }
433  /* option to turn off this element, set abundance to zero */
434  dense.lgElmtOn[nelem] = true;
435  dense.lgElmtSetOff[nelem] = false;
436  /* reset levels to default if they were ever turned off with element off command */
437  if( levels[ipHYDROGEN][nelem] )
438  {
439  iso_sp[ipH_LIKE][nelem].numLevels_max = levels[ipHYDROGEN][nelem];
440  iso_sp[ipHE_LIKE][nelem].numLevels_max = levels[ipHELIUM][nelem];
441  }
442  }
443 
444  else if( p.nMatch("TABL") )
445  {
446 
447  if( !lgElementSet )
448  {
449  fprintf( ioQQQ,
450  " ParseElement did not find an element on the following line:\n" );
451  p.PrintLine( ioQQQ );
453  }
454  /* read in table of position-dependent abundances for a particular element. */
455  abund.lgAbunTabl[nelem] = true;
456 
457  /* general flag saying this option turned on */
458  abund.lgAbTaON = true;
459 
460  /* does the table give depth or radius? keyword DEPTH
461  * says it is depth, default is radius */
462  if( p.nMatch("DEPT") )
463  abund.lgAbTaDepth[nelem] = true;
464  else
465  abund.lgAbTaDepth[nelem] = false;
466 
467  /* make sure not trying to change hydrogen */
468  if( nelem == ipHYDROGEN )
469  {
470  fprintf( ioQQQ, " cannot change abundance of hydrogen.\n" );
471  fprintf( ioQQQ, " Sorry.\n" );
473  }
474 
475  /* read pair giving radius and abundance */
476  p.getline();
477  abund.AbTabRad[0][nelem] = (realnum)p.FFmtRead();
478  abund.AbTabFac[0][nelem] = (realnum)p.FFmtRead();
479 
480  if( p.lgEOL() )
481  {
482  fprintf( ioQQQ, " no pairs entered - cannot interpolate\n" );
484  }
485 
486  /* number of points in the table */
487  abund.nAbunTabl = 2;
488 
489  lgEnd = false;
490  /* LIMTAB is limit to number of points we can store and is
491  * set to 500 in abundances */
492  while( !lgEnd && abund.nAbunTabl < LIMTABD )
493  {
494  p.getline();
495  lgEnd = p.m_lgEOF;
496  if( !lgEnd )
497  {
498  /* convert first 4 char to caps, into chCap */
499  lgEnd = p.hasCommand("END");
500  }
501 
502  /* lgEnd may have been set within above if, if end line encountered*/
503  if( !lgEnd )
504  {
505  abund.AbTabRad[abund.nAbunTabl-1][nelem] =
506  (realnum)p.FFmtRead();
507 
508  abund.AbTabFac[abund.nAbunTabl-1][nelem] =
509  (realnum)p.FFmtRead();
510  abund.nAbunTabl += 1;
511  }
512  }
513 
514  abund.nAbunTabl -= 1;
515 
516  /* now check that radii are in increasing order */
517  for( long i=1; i < abund.nAbunTabl; i++ )
518  {
519  /* the radius values are assumed to be strictly increasing */
520  if( abund.AbTabRad[i][nelem] <= abund.AbTabRad[i-1][nelem] )
521  {
522  fprintf( ioQQQ, "ParseElement: TABLE ELEMENT TABLE radii "
523  "must be in increasing order\n" );
525  }
526  }
527  }
528 
529  else
530  {
531  fprintf( ioQQQ, "ParseElement: ELEMENT command - there must be "
532  "a keyword on this line.\n" );
533  fprintf( ioQQQ, " The keys I know about are TABLE, SCALE, _OFF, "
534  "_ON_, IONIZATION, ABUNDANCE, and ISOTOPES.\n" );
535  fprintf( ioQQQ, " Sorry.\n" );
537  }
538 
539  /* vary option */
540  if( optimize.lgVarOn )
541  {
543  if( p.nMatch("SCAL") )
544  {
545  /* vary scale factor */
546  sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s SCALE %%f LOG",
548 
549  }
550  else if( p.nMatch("ABUN") )
551  {
552  /* vary absolute abundance */
553  sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s ABUNDANCE %%f LOG",
555  }
556  /* pointer to where to write */
558  if( param > 0. )
559  optimize.vparm[0][optimize.nparm] = (realnum)log10(param);
560  else
561  optimize.vparm[0][optimize.nparm] = (realnum)param;
562  optimize.vincr[optimize.nparm] = 0.2f;
564  ++optimize.nparm;
565  }
566  return;
567 }
bool nMatch(const char *chKey) const
Definition: parser.h:150
t_mole_global mole_global
Definition: mole.cpp:7
bool hasCommand(const char *s2)
Definition: parser.cpp:746
bool lgElmtOn
Definition: deuterium.h:21
long int nAbunTabl
Definition: abund.h:235
double FFmtRead(void)
Definition: parser.cpp:472
double exp10(double x)
Definition: cddefines.h:1368
const int ipHE_LIKE
Definition: iso.h:65
vector< bool > lgTreatIsotopes
Definition: mole.h:359
t_input input
Definition: input.cpp:12
long int nvfpnt[LIMPAR]
Definition: optimize.h:198
realnum AbTabFac[LIMTABD][LIMELM]
Definition: abund.h:229
realnum SetIoniz[LIMELM][LIMELM+1]
Definition: dense.h:172
const int NISO
Definition: cddefines.h:311
long int IonHigh[LIMELM+1]
Definition: dense.h:130
long int nRead
Definition: input.h:105
bool lgElmtSetOff[LIMELM]
Definition: dense.h:164
bool isAnyIllegal(void)
Definition: abund.h:145
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:267
bool lgSetIoniz[LIMELM]
Definition: dense.h:167
FILE * ioQQQ
Definition: cddefines.cpp:7
realnum vparm[LIMEXT][LIMPAR]
Definition: optimize.h:192
bool lgTalk
Definition: called.h:12
Definition: parser.h:43
bool lgVarOn
Definition: optimize.h:207
void normAbn()
Definition: abund.h:130
t_dense dense
Definition: global.cpp:15
t_elementnames elementnames
Definition: elementnames.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
bool lgAbnSolar
Definition: abund.h:202
t_abund abund
Definition: abund.cpp:5
t_isotope IsoAbn[LIMELM]
Definition: abund.h:213
bool lgAbTaDepth[LIMELM]
Definition: abund.h:218
long int nparm
Definition: optimize.h:204
int getNiso()
Definition: abund.h:66
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
const realnum BIGFLOAT
Definition: cpu.h:244
#define cdEXIT(FAIL)
Definition: cddefines.h:482
NORETURN void NoNumb(const char *chDesc) const
Definition: parser.cpp:345
long int GetElem(void) const
Definition: parser.cpp:321
const long LIMPAR
Definition: optimize.h:61
t_optimize optimize
Definition: optimize.cpp:6
char chElementNameShort[LIMELM][CHARS_ELEMENT_NAME_SHORT]
Definition: elementnames.h:21
realnum AbTabRad[LIMTABD][LIMELM]
Definition: abund.h:229
realnum vincr[LIMPAR]
Definition: optimize.h:195
#define LIMTABD
Definition: abund.h:226
bool lgElmtOn[LIMELM]
Definition: dense.h:160
long int ipSolar[LIMELM]
Definition: abund.h:238
realnum AbundanceLimit
Definition: dense.h:151
double levels(const genericState &gs)
#define ASSERT(exp)
Definition: cddefines.h:613
realnum getAbn(int Anum)
Definition: abund.h:121
bool getline()
Definition: parser.cpp:273
const int ipH_LIKE
Definition: iso.h:64
bool lgAbTaON
Definition: abund.h:218
const int LIMELM
Definition: cddefines.h:308
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
const int ipHELIUM
Definition: cddefines.h:350
long int npSolar
Definition: abund.h:238
t_deuterium deut
Definition: deuterium.cpp:7
void prtIso(FILE *fp)
Definition: abund.h:158
bool lgEOL(void) const
Definition: parser.h:113
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
long int numLevels_max
Definition: iso.h:524
int PrintLine(FILE *fp) const
Definition: parser.h:206
char chElementName[LIMELM][CHARS_ELEMENT_NAME]
Definition: elementnames.h:17
bool m_lgEOF
Definition: parser.h:55
realnum ScaleElement[LIMELM]
Definition: abund.h:242
int setAbn(int Anum, realnum abn)
Definition: abund.h:114
const int ipHYDROGEN
Definition: cddefines.h:349
void prtIsoPairs(FILE *fp)
Definition: abund.h:169
bool lgAbunTabl[LIMELM]
Definition: abund.h:218
long int nvarxt[LIMPAR]
Definition: optimize.h:198
void ParseElement(Parser &p)
t_called called
Definition: called.cpp:4
realnum solar[LIMELM]
Definition: abund.h:210