cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
abund_starburst.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 /*abund_starburst generate abundance set from Fred Hamann's starburst evolution grid */
4 #include "cddefines.h"
5 #include "optimize.h"
6 #include "input.h"
7 #include "elementnames.h"
8 #include "abund.h"
9 #include "parser.h"
10 
12 {
13  bool lgDebug;
14  long int j;
15  double sqrzed,
16  zHigh,
17  zal,
18  zar,
19  zc,
20  zca,
21  zcl,
22  zco,
23  zed,
24  zed2,
25  zed3,
26  zed4,
27  zedlog,
28  zfe,
29  zh,
30  zhe,
31  zmg,
32  zn,
33  zna,
34  zne,
35  zni,
36  zo,
37  zs,
38  zsi;
39  /* this is largest stored metallicity */
40  static double zLimit = 35.5;
41 
42  DEBUG_ENTRY( "abund_starburst()" );
43 
44  // save default abundances since will use as scale factor
45  for( long i=1; i < LIMELM; i++ )
46  {
47  abund.SolarSave[i] = abund.solar[i];
48  }
49 
50  if( p.nMatch("TRAC") )
51  {
52  lgDebug = true;
53  /* trace keyword did appear
54  * this will not be used, but must trick the sanity test */
55  zHigh = zLimit;
56 
57  /* if trace option set, print header now, and init ZED */
58  fprintf( ioQQQ, " Abundances relative to H, Z\n" );
59 
60  /* this is actual header line */
61  fprintf( ioQQQ, " Z ");
62 
63  /* make line of chemical symbol names */
64  for(j=0; j < 30; j++)
65  {
66  fprintf( ioQQQ, " %2.2s", elementnames.chElementSym[j]);
67  }
68  fprintf( ioQQQ, " \n" );
69 
70  zed = 1e-3;
71  }
72  else
73  {
74  lgDebug = false;
75 
76  /* argument is metallicity */
77  zed = p.getNumberCheckLogLinNegImplLog("metallicity");
78 
79  if( zed <= 0. )
80  {
81  fprintf( ioQQQ, " Z .le.0 not allowed, Z=%10.2e\n",
82  zed );
84  }
85 
86  zHigh = zed;
87  }
88 
89 
90  /* following if loop only if trace option is set
91  * >>chng 97 jun 17, get rid of go to
92  *999 continue
93  * */
94  while( zed <= zHigh )
95  {
96  if( zed < 1e-3 || zed > zLimit )
97  {
98  fprintf( ioQQQ, " The metallicity must be between 1E-3 and%10.2e\n",
99  zLimit );
101  }
102  zed2 = zed*zed;
103  zed3 = zed2*zed;
104  zed4 = zed3*zed;
105  zedlog = log(zed);
106  sqrzed = sqrt(zed);
107  /* The value of each element as a function of ZED=[Z] */
108  zh = 1.081646723 - 0.04600492*zed + 8.6569e-6*zed2 + 1.922e-5*
109  zed3 - 2.0087e-7*zed4;
110 
111  /* helium */
112  zhe = 0.864675891 + 0.044423807*zed + 7.10886e-5*zed2 - 5.3242e-5*
113  zed3 + 5.70194e-7*zed4;
114  abund.solar[1] = (realnum)zhe;
115 
116  /* li, b, bo unchanged */
117  abund.solar[2] = 1.;
118  abund.solar[3] = 1.;
119  abund.solar[4] = 1.;
120 
121  /* carbon */
122  zc = 0.347489799 + 0.016054107*zed - 0.00185855*zed2 + 5.43015e-5*
123  zed3 - 5.3123e-7*zed4;
124  abund.solar[5] = (realnum)zc;
125 
126  /* nitrogen */
127  zn = -0.06549567 + 0.332073984*zed + 0.009146066*zed2 - 0.00054441*
128  zed3 + 6.16363e-6*zed4;
129  if( zn < 0.0 )
130  {
131  zn = 0.05193*zed;
132  }
133  if( zed < 0.3 )
134  {
135  zn = -0.00044731103 + 0.00026453554*zed + 0.52354843*zed2 -
136  0.41156655*zed3 + 0.1290967*zed4;
137  if( zn < 0.0 )
138  {
139  zn = 0.000344828*zed;
140  }
141  }
142  abund.solar[6] = (realnum)zn;
143 
144  /* oxygen */
145  zo = 1.450212747 - 0.05379041*zed + 0.000498919*zed2 + 1.09646e-5*
146  zed3 - 1.8147e-7*zed4;
147  abund.solar[7] = (realnum)zo;
148 
149  /* neon */
150  zne = 1.110139023 + 0.002551998*zed - 2.09516e-7*zed3 - 0.00798157*
151  POW2(zedlog);
152  abund.solar[9] = (realnum)zne;
153 
154  /* fluorine, scale from neon */
155  abund.solar[8] = abund.solar[9];
156 
157  /* sodium */
158  zna = 1.072721387 - 0.02391599*POW2(zedlog) + .068644972*
159  zedlog + 0.017030935/sqrzed;
160  /* this one is slightly negative at very low Z */
161  zna = MAX2(1e-12,zna);
162  abund.solar[10] = (realnum)zna;
163 
164  /* magnesium */
165  zmg = 1.147209646 - 7.9491e-7*POW3(zed) - .00264458*POW2(zedlog) -
166  0.00635552*zedlog;
167  abund.solar[11] = (realnum)zmg;
168 
169  /* aluminium */
170  zal = 1.068116358 - 0.00520227*sqrzed*zedlog - 0.01403851*
171  POW2(zedlog) + 0.066186787*zedlog;
172  /* this one is slightly negative at very low Z */
173  zal = MAX2(1e-12,zal);
174  abund.solar[12] = (realnum)zal;
175 
176  /* silicon */
177  zsi = 1.067372578 + 0.011818743*zed - 0.00107725*zed2 + 3.66056e-5*
178  zed3 - 3.556e-7*zed4;
179  abund.solar[13] = (realnum)zsi;
180 
181  /* phosphorus scaled from silicon */
182  abund.solar[14] = abund.solar[13];
183 
184  /* sulphur */
185  zs = 1.12000;
186  abund.solar[15] = (realnum)zs;
187 
188  /* chlorine */
189  zcl = 1.10000;
190  abund.solar[16] = (realnum)zcl;
191 
192  /* argon */
193  zar = 1.091067724 + 2.51124e-6*zed3 - 0.0039589*sqrzed*zedlog +
194  0.015686715*zedlog;
195  abund.solar[17] = (realnum)zar;
196 
197  /* potassium scaled from silicon */
198  abund.solar[18] = abund.solar[13];
199 
200  /* calcium */
201  zca = 1.077553875 - 0.00888806*zed + 0.001479866*zed2 - 6.5689e-5*
202  zed3 + 1.16935e-6*zed4;
203  abund.solar[19] = (realnum)zca;
204 
205  /* iron */
206  zfe = 0.223713045 + 0.001400746*zed + 0.000624734*zed2 - 3.5629e-5*
207  zed3 + 8.13184e-7*zed4;
208  abund.solar[25] = (realnum)zfe;
209 
210  /* scandium, scaled from iron */
211  abund.solar[20] = abund.solar[25];
212 
213  /* titanium, scaled from iron */
214  abund.solar[21] = abund.solar[25];
215 
216  /* vanadium scaled from iron */
217  abund.solar[22] = abund.solar[25];
218 
219  /* chromium scaled from iron */
220  abund.solar[23] = abund.solar[25];
221 
222  /* manganese scaled from iron */
223  abund.solar[24] = abund.solar[25];
224 
225  /* cobalt */
226  zco = zfe;
227  abund.solar[26] = (realnum)zco;
228 
229  /* nickel */
230  zni = zfe;
231  abund.solar[27] = (realnum)zni;
232 
233  /* copper scaled from iron */
234  abund.solar[28] = abund.solar[25];
235 
236  /* zinc scaled from iron */
237  abund.solar[29] = abund.solar[25];
238 
239  /* rescale to true abundances */
240  abund.solar[0] = 1.;
242  zh);
243 
244  for( long i=2; i < LIMELM; i++ )
245  {
247  zed/zh);
248  }
249 
250  if( lgDebug )
251  {
252  fprintf( ioQQQ, "%10.2e", zed );
253  for( long i=0; i < LIMELM; i++ )
254  {
255  fprintf( ioQQQ, "%6.2f", log10(abund.solar[i]) );
256  }
257  fprintf( ioQQQ, "\n" );
258 
259  if( fp_equal( zed, zLimit ) )
260  {
262  }
263  }
264 
265  /* this trick makes sure last one is upper limit */
266  if( zed < zLimit )
267  {
268  zed = MIN2(zed*1.5,zLimit);
269  }
270  else
271  {
272  zed *= 1.5;
273  }
274  }
275 
276  /* vary option */
277  if( optimize.lgVarOn )
278  {
279  /* this is number of parameters to feed onto the input line */
281  strcpy( optimize.chVarFmt[optimize.nparm], "ABUNDANCES STARBURST %f LOG" );
282  /* following are upper and lower limits to metallicity range */
283  optimize.varang[optimize.nparm][0] = (realnum)log10(1e-3);
284  optimize.varang[optimize.nparm][1] = (realnum)log10(zLimit);
285  /* pointer to where to write */
287  /* log of metallicity will be argument */
288  optimize.vparm[0][optimize.nparm] = (realnum)log10(zed);
289  optimize.vincr[optimize.nparm] = 0.2f;
290  ++optimize.nparm;
291  }
292  return;
293 }
bool nMatch(const char *chKey) const
Definition: parser.h:150
t_input input
Definition: input.cpp:12
double getNumberCheckLogLinNegImplLog(const char *chDesc)
Definition: parser.cpp:410
long int nvfpnt[LIMPAR]
Definition: optimize.h:198
realnum varang[LIMPAR][2]
Definition: optimize.h:201
long int nRead
Definition: input.h:105
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:267
FILE * ioQQQ
Definition: cddefines.cpp:7
realnum vparm[LIMEXT][LIMPAR]
Definition: optimize.h:192
realnum SolarSave[LIMELM]
Definition: abund.h:200
#define MIN2(a, b)
Definition: cddefines.h:803
Definition: parser.h:43
bool lgVarOn
Definition: optimize.h:207
t_elementnames elementnames
Definition: elementnames.cpp:5
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:854
t_abund abund
Definition: abund.cpp:5
#define POW2
Definition: cddefines.h:979
long int nparm
Definition: optimize.h:204
float realnum
Definition: cddefines.h:124
#define EXIT_FAILURE
Definition: cddefines.h:168
#define cdEXIT(FAIL)
Definition: cddefines.h:482
t_optimize optimize
Definition: optimize.cpp:6
realnum vincr[LIMPAR]
Definition: optimize.h:195
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
const int LIMELM
Definition: cddefines.h:308
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
#define MAX2(a, b)
Definition: cddefines.h:824
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
void abund_starburst(Parser &p)
#define POW3
Definition: cddefines.h:986
long int nvarxt[LIMPAR]
Definition: optimize.h:198
realnum solar[LIMELM]
Definition: abund.h:210