cloudy  trunk
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 
11 void abund_starburst(Parser &p)
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 );
77  cdEXIT(EXIT_FAILURE);
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 );
94  cdEXIT(EXIT_FAILURE);
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  {
255  cdEXIT(EXIT_FAILURE);
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 }
Parser::nMatch
bool nMatch(const char *chKey) const
Definition: parser.h:135
t_optimize::vincr
realnum vincr[LIMPAR]
Definition: optimize.h:191
t_optimize::nparm
long int nparm
Definition: optimize.h:201
elementnames.h
t_input::nRead
long int nRead
Definition: input.h:49
ioQQQ
FILE * ioQQQ
Definition: cddefines.cpp:7
elementnames
t_elementnames elementnames
Definition: elementnames.cpp:5
realnum
float realnum
Definition: cddefines.h:103
abund.h
t_optimize::lgVarOn
bool lgVarOn
Definition: optimize.h:203
Parser::getNumberCheckLogLinNegImplLog
double getNumberCheckLogLinNegImplLog(const char *chDesc)
Definition: parser.cpp:291
input
t_input input
Definition: input.cpp:12
POW2
#define POW2
Definition: cddefines.h:929
MIN2
#define MIN2
Definition: cddefines.h:761
optimize
t_optimize optimize
Definition: optimize.cpp:5
t_optimize::vparm
realnum vparm[LIMEXT][LIMPAR]
Definition: optimize.h:188
t_optimize::chVarFmt
char chVarFmt[LIMPAR][FILENAME_PATH_LENGTH_2]
Definition: optimize.h:263
EXIT_FAILURE
#define EXIT_FAILURE
Definition: cddefines.h:140
Parser
Definition: parser.h:31
cddefines.h
POW3
#define POW3
Definition: cddefines.h:936
optimize.h
t_abund::solar
realnum solar[LIMELM]
Definition: abund.h:65
t_elementnames::chElementSym
char chElementSym[LIMELM][CHARS_ELEMENT_SYM]
Definition: elementnames.h:25
t_optimize::nvarxt
long int nvarxt[LIMPAR]
Definition: optimize.h:194
MAX2
#define MAX2
Definition: cddefines.h:782
LIMELM
const int LIMELM
Definition: cddefines.h:258
cdEXIT
#define cdEXIT(FAIL)
Definition: cddefines.h:434
abund
t_abund abund
Definition: abund.cpp:5
t_optimize::varang
realnum varang[LIMPAR][2]
Definition: optimize.h:198
fp_equal
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition: cddefines.h:812
parser.h
t_optimize::nvfpnt
long int nvfpnt[LIMPAR]
Definition: optimize.h:195
abund_starburst
void abund_starburst(Parser &p)
Definition: abund_starburst.cpp:11
t_abund::SolarSave
realnum SolarSave[LIMELM]
Definition: abund.h:46
input.h
DEBUG_ENTRY
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:684
Generated by   doxygen 1.8.17