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