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