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

Generated for cloudy by doxygen 1.8.1.1