cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
abundances.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 /*AbundancesPrt print all abundances, both gas phase and grains */
4 /*AbundancesSet sets initial abundances after parameters are entered by reading input */
5 /*AbundancesTable interpolate on table of points to do 'element table' command, */
6 /*PrtElem print chemical composition at start of calculation */
7 #include "cddefines.h"
8 #include "phycon.h"
9 #include "called.h"
10 #include "stopcalc.h"
11 #include "thermal.h"
12 #include "trace.h"
13 #include "elementnames.h"
14 #include "dense.h"
15 #include "radius.h"
16 #include "grainvar.h"
17 #include "abund.h"
18 #include "deuterium.h"
19 #include "physconst.h"
20 
21 /*PrtElem print chemical composition at start of calculation */
22 STATIC void PrtElem(
23  /* the job to do, the options are "init", "fill", "flus" */
24  const char *chJob,
25  /* label for the element */
26  const char *chLabl,
27  /* its abundance */
28  double abund_prt);
29 
30 /*AbundancesPrt print all abundances, both gas phase and grains */
31 void AbundancesPrt( void )
32 {
33  long int i;
34  double GrainNumRelHydrSilicate ,
35  GrainNumRelHydrCarbonaceous ,
36  GrainNumRelHydr_PAH,
37  GrainMassRelHydrSilicate,
38  GrainMassRelHydrCarbonaceous,
39  GrainMassRelHydr_PAH;
40 
41  DEBUG_ENTRY( "AbundancesPrt()" );
42 
43  /* this is main loop to print abundances of each element */
44  if( called.lgTalk )
45  {
46  PrtElem("initG"," ",0.);/* initialize print routine for gas*/
47  for( i=0; i < LIMELM; i++ )
48  {
49  if( dense.lgElmtOn[i] )
50  {
51  /* fill in print buffer with abundances */
52  PrtElem("fill",(char*)elementnames.chElementSym[i],
53  abund.solar[i]);
54  }
55  }
56 
57  /* flush the print buffer */
58  PrtElem("flus"," ",0.);
59  /* final carriage return */
60  fprintf( ioQQQ, " \n" );
61 
62  /* now grains if present */
63  if( gv.lgDustOn() )
64  {
65  /* we will first print the total abundances of each element locked up in grains */
66  /* initialize print routine for dust*/
67  PrtElem("initD"," ",0.);
68  for( i=0; i < LIMELM; i++ )
69  {
70  if( gv.elmSumAbund[i]>SMALLFLOAT )
71  {
72  /* fill in print buffer with abundances */
73  PrtElem("fill",(char*)elementnames.chElementSym[i],
75  }
76  }
77  /* flush the print buffer */
78  PrtElem("flus"," ",0.);
79  /* final carriage return */
80  fprintf( ioQQQ, " \n" );
81 
82  /* this is used to store grain number density per hydrogen */
83  GrainNumRelHydrSilicate = 0.;
84  GrainNumRelHydrCarbonaceous = 0;
85  GrainNumRelHydr_PAH = 0.;
86  GrainMassRelHydrSilicate = 0.;
87  GrainMassRelHydrCarbonaceous = 0;
88  GrainMassRelHydr_PAH = 0.;
89 
90  for( size_t nd=0; nd < gv.bin.size(); nd++ )
91  {
92 
93  /* number density of grains per hydrogen, the ratio
94  * gv.bin[nd]->IntVol/gv.bin[nd]->AvVol is the number of grain particles
95  * per H at standard grain abundance*/
96  realnum DensityNumberPerHydrogen =
97  (gv.bin[nd]->IntVol/gv.bin[nd]->AvVol)*gv.bin[nd]->dstAbund /
98  gv.bin[nd]->GrnDpth;
99  /* mass of grains per hydrogen */
100  realnum DensityMassPerHydrogen =
101  gv.bin[nd]->IntVol*gv.bin[nd]->dustp[0]*gv.bin[nd]->dstAbund/
102  (realnum)ATOMIC_MASS_UNIT / gv.bin[nd]->GrnDpth;
103 
104  /* >>chng 06 mar 05, fix expression for calculating grain number density, PvH */
105  if( gv.bin[nd]->matType == MAT_CAR || gv.bin[nd]->matType == MAT_CAR2 ||
106  gv.bin[nd]->matType == MAT_SIC )
107  {
108  /* carbonaceous grains */
109  GrainNumRelHydrCarbonaceous += DensityNumberPerHydrogen;
110  GrainMassRelHydrCarbonaceous += DensityMassPerHydrogen;
111  }
112  else if( gv.bin[nd]->matType == MAT_SIL || gv.bin[nd]->matType == MAT_SIL2 )
113  {
114  /* silicate grains */
115  GrainNumRelHydrSilicate += DensityNumberPerHydrogen;
116  GrainMassRelHydrSilicate += DensityMassPerHydrogen;
117  }
118  else if( gv.bin[nd]->matType == MAT_PAH || gv.bin[nd]->matType == MAT_PAH2 )
119  {
120  /* PAHs - full abundance - remove possible factor accounting for
121  * variation of abundances with physical conditions - this will
122  * be the PAH abundance with scale factor of unity */
123  GrainNumRelHydr_PAH += DensityNumberPerHydrogen;
124  GrainMassRelHydr_PAH += DensityMassPerHydrogen;
125  }
126  else
127  TotalInsanity();
128  }
129 
130  /* now print total number of grains of each type */
131  fprintf(ioQQQ," Number of grains per hydrogen (scale=1) Mass of grains per hydrogen (scale=1)\n");
132  fprintf(ioQQQ," Carbonaceous: %.3f Silicate: %.3f PAH: %.3f Carbonaceous: %.3f Silicate: %.3f PAH: %.3f\n\n" ,
133  log10( MAX2( 1e-30, GrainNumRelHydrCarbonaceous ) ) ,
134  log10( MAX2( 1e-30, GrainNumRelHydrSilicate ) ) ,
135  log10( MAX2( 1e-30, GrainNumRelHydr_PAH ) ) ,
136  log10( MAX2( 1e-30, GrainMassRelHydrCarbonaceous ) ) ,
137  log10( MAX2( 1e-30, GrainMassRelHydrSilicate ) ) ,
138  log10( MAX2( 1e-30, GrainMassRelHydr_PAH ) ) );
139  }
140  }
141  return;
142 }
143 
144 /*AbundancesSet print all abundances, both gas phase and grains */
145 void AbundancesSet(void)
146 {
147  long int i,
148  nelem;
149  double fac;
150  static bool lgFirstCall=true;
151  static bool lgElOnOff[LIMELM];
152 
153  DEBUG_ENTRY( "AbundancesSet()" );
154 
155  /* if this is the first call to this routine in this core load,
156  * save the state of the lgElmOn array, so that it is possible
157  * to turn off elements in later models, but not turn on an
158  * element that was initially turned off. This is necessary since
159  * the Create... routines that create space for elements will
160  * not be revisited in later models. You can turn off an initially
161  * enabled element, but not turn a disabled one on. */
162 
163  if( lgFirstCall )
164  {
165  /* first call - save the initial state of the lgElmtOn vector */
166  for( i=0; i<LIMELM; ++i )
167  {
168  lgElOnOff[i] = dense.lgElmtOn[i];
169  }
170  }
171  lgFirstCall = false;
172 
173  /* make sure that initially false elements remain off, while letting
174  * enabled elements be turned off */
175  for( i=ipHYDROGEN; i<LIMELM; ++i )
176  {
177  dense.lgElmtOn[i] = lgElOnOff[i] && dense.lgElmtOn[i];
178  }
179 
180  /* rescale so that abundances are H=1 */
181  for( i=ipHELIUM; i < LIMELM; i++ )
182  {
183  abund.solar[i] /= abund.solar[0];
184  }
185  abund.solar[ipHYDROGEN] = 1.;
186 
187  /* set current abundances to "solar" times metals scale factor
188  * and grain depletion factor */
190 
191  /* option for density or abundance variations, this flag is true by default,
192  * set in zero, but set false if variations are enabled AND these
193  * are not density variations, but rather abundances */
194  if( dense.lgDenFlucOn )
195  {
196  /* usual case - either density fluctuations or none at all */
197  fac = 1.;
198  }
199  else
200  {
201  /* abundance fluctuations enabled, set initial value */
202  fac = dense.cfirst*cos(dense.flcPhase) + dense.csecnd;
203  }
204 
205  for( i=ipLITHIUM; i < LIMELM; i++ )
206  {
208  abund.ScaleElement[i]*fac);
209  }
210 
211  /* now fix abundance of any element with element table set */
212  if( abund.lgAbTaON )
213  {
214  for( nelem=ipHELIUM; nelem < LIMELM; ++nelem )
215  {
216  if( abund.lgAbunTabl[nelem] )
217  {
219  radius.depth,nelem+1));
220  }
221  }
222  }
223 
224  /* dense.gas_phase[nelem] contains total abundance of element */
225  /* the density of hydrogen itself has already been set at this point -
226  * it is set when commands parsed, most likely by the hden command -
227  * set all heavier elements */
228  /* if abund.solar[ipHYDROGEN] == 1, consistency doesn't hurt that much */
229  for( nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
230  {
231  /* this implements the element off limit xxx command, where
232  * xxx is the limit to the smallest n(A)/n(H) that will remain on */
233  if( abund.solar[nelem] < dense.AbundanceLimit )
234  dense.lgElmtOn[nelem] = false;
235 
236  if( dense.lgElmtOn[nelem] )
237  {
239  if( dense.gas_phase[nelem] <= 0. )
240  {
241  fprintf( ioQQQ, " Abundances must be greater than zero. "
242  "Check entered abundance for element%3ld = %2.2s, value=%.2e\n",
243  nelem, elementnames.chElementSym[nelem] , dense.gas_phase[nelem]);
245  }
246  else if( dense.gas_phase[nelem] < SMALLFLOAT )
247  {
248  fprintf(ioQQQ," Abundance for %s is %.2e, less than lower "
249  "limit of %.3e, so turning element off.\n",
250  elementnames.chElementSym[nelem],
251  dense.gas_phase[nelem],
252  SMALLFLOAT );
253  dense.lgElmtOn[nelem] = false;
254  }
255  else if( dense.gas_phase[nelem] > MAX_DENSITY )
256  {
257  fprintf(ioQQQ," Abundance for %s is %.2e. This version of Cloudy does not "
258  "permit densities greater than %e cm-3.\n",
259  elementnames.chElementSym[nelem],
260  dense.gas_phase[nelem],
261  MAX_DENSITY );
263  }
264  }
265  if( !dense.lgElmtOn[nelem] )
266  {
267  /* >>chng 04 apr 20, set to zero if element is off */
268  dense.SetGasPhaseDensity( nelem, 0. );
269  }
270 
271  /* Set all neutral ions to maintain invariant */
272  dense.xIonDense[nelem][0] = dense.gas_phase[nelem];
273  for( long int ion=1; ion < LIMELM+1; ion++ )
274  {
275  dense.xIonDense[nelem][ion] = 0.;
276  }
277 
278  if( deut.lgElmtOn )
280  }
281 
282  SumDensities();
283 
284  /* if stop temp set below default then we are going into cold and possibly
285  * molecular gas - check some parameters in this case */
287  /* thermal.ConstTemp def is zero, set pos when used */
289  {
290 
291  /* print warning if temperature set below default but C > O */
293  {
294  fprintf( ioQQQ, "\n >>> \n"
295  " >>> The simulation is going into possibly molecular gas but the carbon/oxygen abundance ratio is greater than unity.\n" );
296  fprintf( ioQQQ, " >>> Standard interstellar chemistry networks are designed for environments with C/O < 1.\n" );
297  fprintf( ioQQQ, " >>> The chemistry network may (or may not) collapse deep in molecular regions where CO is fully formed.\n" );
298  fprintf( ioQQQ, " >>> \n\n\n\n\n" );
299  }
300  }
301 
302  if( trace.lgTrace )
303  {
304  realnum sumx , sumy , sumz = 0.;
305 
308 
309  fprintf( ioQQQ, "\n AbundancesSet sets following densities (cm^-3); \n" );
310  for( i=0; i<3; i++ )
311  {
312  for( nelem=i*10; nelem < i*10+10; nelem++ )
313  {
314  fprintf( ioQQQ, " %2.2s", elementnames.chElementSym[nelem] );
315  PrintE82( ioQQQ, dense.gas_phase[nelem] );
316  if( nelem>ipHELIUM )
317  sumz += dense.gas_phase[nelem]*dense.AtomicWeight[nelem];
318  }
319  fprintf( ioQQQ, " \n" );
320  }
321  fprintf( ioQQQ, "\n AbundancesSet sets following abundances rel to H; \n" );
322  for( i=0; i<3; i++ )
323  {
324  for( nelem=i*10; nelem < i*10+10; nelem++ )
325  {
326  fprintf( ioQQQ, " %2.2s", elementnames.chElementSym[nelem] );
328  }
329  fprintf( ioQQQ, " \n" );
330  }
331  fprintf( ioQQQ, " \n" );
332  fprintf(ioQQQ," Gas-phase mass fractions, X:%.3e Y:%.3e Z:%.3e\n\n",
333  sumx/SDIV(sumx+sumy+sumz) ,
334  sumy/SDIV(sumx+sumy+sumz) ,
335  sumz/SDIV(sumx+sumy+sumz) );
336  }
337  return;
338 }
339 
340 /* this is number of elements across one line */
341 #define NELEM1LINE 9
342 
343 /*PrtElem print chemical composition at start of calculation */
345  /* the job to do, the options are "init", "fill", "flus" */
346  const char *chJob,
347  /* label for the element */
348  const char *chLabl,
349  /* its abundance */
350  double abund_prt)
351 {
352  static char chAllLabels[NELEM1LINE][14];/* buffer where elements will be stored*/
353  long int i,
354  noffset;
355  static long int nelem; /* counter for number of elements read in*/
356 
357  DEBUG_ENTRY( "PrtElem()" );
358 
359  if( strcmp(chJob,"initG") == 0 )
360  {
361  /* gas phase abundances */
362  nelem = 0;
363  fprintf( ioQQQ,
364  " Gas Phase Chemical Composition\n" );
365  }
366  else if( strcmp(chJob,"initD") == 0 )
367  {
368  /* abundances in grains */
369  nelem = 0;
370  fprintf( ioQQQ,
371  " Grain Chemical Composition\n" );
372  }
373 
374  else if( strcmp(chJob,"fill") == 0 )
375  {
376  /* print log of abundance to avoid exponential output */
377  abund_prt = log10( abund_prt );
378  /* stuff in labels and abundances */
379  sprintf( chAllLabels[nelem], " %2.2s:%8.4f", chLabl, abund_prt );
380  if( nelem == NELEM1LINE-1 )
381  {
382  /* we hit as many as it will hold - print it out and reset*/
383  fprintf( ioQQQ, " " );
384  for( i=0; i < NELEM1LINE; i++ )
385  {
386  fprintf( ioQQQ, "%13.13s", chAllLabels[i] );
387  }
388  fprintf( ioQQQ, "\n" );
389  /* reset counter to zero */
390  nelem = 0;
391  }
392  else
393  {
394  /* just increment */
395  ++nelem;
396  }
397  }
398 
399 # if 0
400  /* Do this if you want to know about PAH number abundance */
401  else if( strcmp(chJob,"fillp") == 0 )
402  {
403  /* print log of abundance to avoid exponential output */
404  abund_prt = log10( abund_prt );
405 
406  /* stuff in labels and abundances */
407  sprintf( chAllLabels[nelem], " %2.2s:%8.4f", chLabl, abund_prt );
408  if( nelem == NELEM1LINE-1 )
409  {
410  /* we hit as many as it will hold - print it out and reset*/
411  fprintf( ioQQQ, " " );
412  for( i=0; i < NELEM1LINE; i++ )
413  {
414  fprintf( ioQQQ, "%13.13s", chAllLabels[i] );
415  }
416  fprintf( ioQQQ, "\n" );
417  /* reset counter to zero */
418  nelem = 0;
419  }
420  else
421  {
422  /* just increment */
423  ++nelem;
424  }
425  }
426 # endif
427 
428  else if( strcmp(chJob,"flus") == 0 )
429  {
430  /* flush the stack */
431  i = NELEM1LINE - (nelem - 2);
432  noffset = i/2-1;
433  /* make format pretty */
434  fprintf( ioQQQ, " " );
435 
436  for(i=0; i < noffset; i++)
437  {
438  /* skip out this many fields */
439  fprintf( ioQQQ, " " );
440  }
441 
442  /* if nelem is even we need to space out another 8 */
443  if( !(nelem%2) && nelem > 0)
444  fprintf( ioQQQ," ");
445 
446  for( i=0; i < nelem; i++ )
447  {
448  fprintf( ioQQQ, "%13.13s", chAllLabels[i] );
449  }
450 
451  fprintf( ioQQQ, "\n" );
452  }
453  else
454  {
455  fprintf( ioQQQ, " PrtElem does not understand job=%4.4s\n",
456  chJob );
458  }
459  return;
460 }
461 
462 
463 /*AbundancesTable interpolate on table of points to do 'element table' command, */
464 double AbundancesTable(double r0,
465  double depth,
466  long int iel)
467 {
468  bool lgHit;
469  long int j;
470  double frac,
471  tababun_v,
472  x;
473 
474  DEBUG_ENTRY( "AbundancesTable()" );
475  /* interpolate on table of points to do 'element table' command, based
476  * on code by K Volk, each line is log radius and abundance. */
477 
478  /* interpolate on radius or depth? */
479  if( abund.lgAbTaDepth[iel-1] )
480  {
481  /* depth key appeared = we want depth */
482  x = log10(depth);
483  }
484  else
485  {
486  /* use radius */
487  x = log10(r0);
488  }
489 
490  /* this will be reset below, but is here as a safety check */
491  tababun_v = -DBL_MAX;
492 
493  if( x < abund.AbTabRad[0][iel-1] || x >= abund.AbTabRad[abund.nAbunTabl-1][iel-1] )
494  {
495  fprintf( ioQQQ, " requested radius outside range of AbundancesTable\n" );
496  fprintf( ioQQQ, " radius was%10.2e min, max=%10.2e%10.2e\n",
497  x, abund.AbTabRad[0][iel-1], abund.AbTabRad[abund.nAbunTabl-1][iel-1] );
499  }
500 
501  else
502  {
503  lgHit = false;
504  j = 1;
505 
506  while( !lgHit && j <= abund.nAbunTabl - 1 )
507  {
508  if( abund.AbTabRad[j-1][iel-1] <= (realnum)x &&
509  abund.AbTabRad[j][iel-1] > (realnum)x )
510  {
511  frac = (x - abund.AbTabRad[j-1][iel-1])/(abund.AbTabRad[j][iel-1] -
512  abund.AbTabRad[j-1][iel-1]);
513  tababun_v = abund.AbTabFac[j-1][iel-1] + frac*
514  (abund.AbTabFac[j][iel-1] - abund.AbTabFac[j-1][iel-1]);
515  lgHit = true;
516  }
517  ++j;
518  }
519 
520  if( !lgHit )
521  {
522  fprintf( ioQQQ, " radius outran dlaw table scale, requested=%6.2f largest=%6.2f\n",
523  x, abund.AbTabRad[abund.nAbunTabl-1][iel-1] );
525  }
526  }
527 
528  /* got it, now return value, not log of density */
529  tababun_v = exp10(tababun_v);
530  return( tababun_v );
531 }
532 
double Radius
Definition: radius.h:31
double depth
Definition: radius.h:31
t_thermal thermal
Definition: thermal.cpp:6
void SetGasPhaseDensity(const long nelem, const realnum density)
Definition: dense.cpp:106
bool lgElmtOn
Definition: deuterium.h:21
long int nAbunTabl
Definition: abund.h:235
realnum depset[LIMELM]
Definition: abund.h:245
realnum flcPhase
Definition: dense.h:265
realnum csecnd
Definition: dense.h:264
#define NELEM1LINE
Definition: abundances.cpp:341
double exp10(double x)
Definition: cddefines.h:1368
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
realnum AbTabFac[LIMTABD][LIMELM]
Definition: abund.h:229
const realnum SMALLFLOAT
Definition: cpu.h:246
const int ipOXYGEN
Definition: cddefines.h:356
void AbundancesSet(void)
Definition: abundances.cpp:145
t_StopCalc StopCalc
Definition: stopcalc.cpp:7
STATIC void PrtElem(const char *chJob, const char *chLabl, double abund_prt)
Definition: abundances.cpp:344
t_phycon phycon
Definition: phycon.cpp:6
realnum elmSumAbund[LIMELM]
Definition: grainvar.h:510
FILE * ioQQQ
Definition: cddefines.cpp:7
bool lgTalk
Definition: called.h:12
void InitDeuteriumIonization()
Definition: deuterium.cpp:29
t_dense dense
Definition: global.cpp:15
const double MAX_DENSITY
Definition: cddefines.h:319
t_elementnames elementnames
Definition: elementnames.cpp:5
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
t_trace trace
Definition: trace.cpp:5
t_abund abund
Definition: abund.cpp:5
const double TEMP_STOP_DEFAULT
Definition: phycon.h:119
realnum ScaleMetals
Definition: abund.h:254
void AbundancesPrt(void)
Definition: abundances.cpp:31
bool lgAbTaDepth[LIMELM]
Definition: abund.h:218
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
#define cdEXIT(FAIL)
Definition: cddefines.h:482
double AbundancesTable(double r0, double depth, long int iel)
Definition: abundances.cpp:464
realnum AbTabRad[LIMTABD][LIMELM]
Definition: abund.h:229
t_radius radius
Definition: radius.cpp:5
realnum TempLoStopZone
Definition: stopcalc.h:42
bool lgElmtOn[LIMELM]
Definition: dense.h:160
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
realnum gas_phase[LIMELM]
Definition: dense.h:76
realnum AbundanceLimit
Definition: dense.h:151
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
realnum ConstTemp
Definition: thermal.h:56
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
vector< GrainBin * > bin
Definition: grainvar.h:583
t_deuterium deut
Definition: deuterium.cpp:7
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
bool lgDenFlucOn
Definition: dense.h:255
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
const int ipCARBON
Definition: cddefines.h:354
void SumDensities(void)
Definition: dense.cpp:235
GrainVar gv
Definition: grainvar.cpp:5
realnum ScaleElement[LIMELM]
Definition: abund.h:242
void PrintE82(FILE *, double)
Definition: service.cpp:824
double frac(double d)
realnum cfirst
Definition: dense.h:263
bool lgDustOn() const
Definition: grainvar.h:475
const int ipHYDROGEN
Definition: cddefines.h:349
const int ipLITHIUM
Definition: cddefines.h:351
bool lgAbunTabl[LIMELM]
Definition: abund.h:218
t_called called
Definition: called.cpp:4
realnum solar[LIMELM]
Definition: abund.h:210