root/branches/newmole/source/abund_starburst.cpp

Revision 1780, 7.0 kB (checked in by rjrw, 11 months ago)

Merged from trunk r1738:1779

  • Property svn:eol-style set to native
Line 
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
10void 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.;
269                abund.solar[1] = (realnum)(abund.solar[1]*abund.SolarSave[1]/
270                  zh);
271
272                for( i=2; i < LIMELM; i++ )
273                {
274                        abund.solar[i] = (realnum)(abund.solar[i]*abund.SolarSave[i]*
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 */
308                optimize.nvarxt[optimize.nparm] = 1;
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 */
314                optimize.nvfpnt[optimize.nparm] = input.nRead;
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}
Note: See TracBrowser for help on using the browser.