00001
00002
00003
00004 #include "cddefines.h"
00005 #include "optimize.h"
00006 #include "input.h"
00007 #include "elementnames.h"
00008 #include "abund.h"
00009
00010 void abund_starburst(char *chCard)
00011 {
00012 bool lgDebug,
00013 lgEOL;
00014 long int i,
00015 j;
00016 double sqrzed,
00017 zHigh,
00018 zal,
00019 zar,
00020 zc,
00021 zca,
00022 zcl,
00023 zco,
00024 zed,
00025 zed2,
00026 zed3,
00027 zed4,
00028 zedlog,
00029 zfe,
00030 zh,
00031 zhe,
00032 zmg,
00033 zn,
00034 zna,
00035 zne,
00036 zni,
00037 zo,
00038 zs,
00039 zsi;
00040
00041 static double zLimit = 35.5;
00042
00043 DEBUG_ENTRY( "abund_starburst()" );
00044
00045
00046 if( nMatch("TRAC",chCard) )
00047 {
00048 lgDebug = true;
00049 }
00050 else
00051 {
00052 lgDebug = false;
00053 }
00054
00055 i = 5;
00056
00057 zed = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
00058 if( lgDebug )
00059 {
00060
00061
00062 zHigh = zLimit;
00063
00064
00065 fprintf( ioQQQ, " Abundances relative to H, Z\n" );
00066
00067
00068 fprintf( ioQQQ, " Z ");
00069
00070
00071 for(j=0; j < 30; j++)
00072 {
00073 fprintf( ioQQQ, " %2.2s", elementnames.chElementSym[j]);
00074 }
00075 fprintf( ioQQQ, " \n" );
00076
00077 zed = 1e-3;
00078 }
00079 else if( lgEOL && !lgDebug )
00080 {
00081
00082 fprintf( ioQQQ, " The metallicity must appear on this line.\n" );
00083 puts( "[Stop in abund_starburst]" );
00084 cdEXIT(EXIT_FAILURE);
00085 }
00086 else
00087 {
00088 zHigh = zed;
00089 }
00090
00091
00092
00093 if( nMatch(" LOG",chCard) )
00094 {
00095 zed = pow(10.,zed);
00096 zHigh = zed;
00097 }
00098
00099 else if( nMatch("LINE",chCard) )
00100 {
00101 if( zed <= 0. )
00102 {
00103 fprintf( ioQQQ, " Z .le.0 not allowed, Z=%10.2e\n",
00104 zed );
00105 puts( "[Stop in abund_starburst]" );
00106 cdEXIT(EXIT_FAILURE);
00107 }
00108 }
00109
00110 else
00111 {
00112
00113 if( zed <= 0. )
00114 {
00115 zed = pow(10.,zed);
00116 zHigh = zed;
00117 }
00118 }
00119
00120
00121
00122
00123
00124 while( zed <= zHigh )
00125 {
00126 if( zed < 1e-3 || zed > zLimit )
00127 {
00128 fprintf( ioQQQ, " The metallicity must be between 1E-3 and%10.2e\n",
00129 zLimit );
00130 puts( "[Stop in abund_starburst]" );
00131 cdEXIT(EXIT_FAILURE);
00132 }
00133 zed2 = zed*zed;
00134 zed3 = zed2*zed;
00135 zed4 = zed3*zed;
00136 zedlog = log(zed);
00137 sqrzed = sqrt(zed);
00138
00139 zh = 1.081646723 - 0.04600492*zed + 8.6569e-6*zed2 + 1.922e-5*
00140 zed3 - 2.0087e-7*zed4;
00141
00142
00143 zhe = 0.864675891 + 0.044423807*zed + 7.10886e-5*zed2 - 5.3242e-5*
00144 zed3 + 5.70194e-7*zed4;
00145 abund.solar[1] = (float)zhe;
00146
00147
00148 abund.solar[2] = 1.;
00149 abund.solar[3] = 1.;
00150 abund.solar[4] = 1.;
00151
00152
00153 zc = 0.347489799 + 0.016054107*zed - 0.00185855*zed2 + 5.43015e-5*
00154 zed3 - 5.3123e-7*zed4;
00155 abund.solar[5] = (float)zc;
00156
00157
00158 zn = -0.06549567 + 0.332073984*zed + 0.009146066*zed2 - 0.00054441*
00159 zed3 + 6.16363e-6*zed4;
00160 if( zn < 0.0 )
00161 {
00162 zn = 0.05193*zed;
00163 }
00164 if( zed < 0.3 )
00165 {
00166 zn = -0.00044731103 + 0.00026453554*zed + 0.52354843*zed2 -
00167 0.41156655*zed3 + 0.1290967*zed4;
00168 if( zn < 0.0 )
00169 {
00170 zn = 0.000344828*zed;
00171 }
00172 }
00173 abund.solar[6] = (float)zn;
00174
00175
00176 zo = 1.450212747 - 0.05379041*zed + 0.000498919*zed2 + 1.09646e-5*
00177 zed3 - 1.8147e-7*zed4;
00178 abund.solar[7] = (float)zo;
00179
00180
00181 zne = 1.110139023 + 0.002551998*zed - 2.09516e-7*zed3 - 0.00798157*
00182 POW2(zedlog);
00183 abund.solar[9] = (float)zne;
00184
00185
00186 abund.solar[8] = abund.solar[9];
00187
00188
00189 zna = 1.072721387 - 0.02391599*POW2(zedlog) + .068644972*
00190 zedlog + 0.017030935/sqrzed;
00191
00192 zna = MAX2(1e-12,zna);
00193 abund.solar[10] = (float)zna;
00194
00195
00196 zmg = 1.147209646 - 7.9491e-7*POW3(zed) - .00264458*POW2(zedlog) -
00197 0.00635552*zedlog;
00198 abund.solar[11] = (float)zmg;
00199
00200
00201 zal = 1.068116358 - 0.00520227*sqrzed*zedlog - 0.01403851*
00202 POW2(zedlog) + 0.066186787*zedlog;
00203
00204 zal = MAX2(1e-12,zal);
00205 abund.solar[12] = (float)zal;
00206
00207
00208 zsi = 1.067372578 + 0.011818743*zed - 0.00107725*zed2 + 3.66056e-5*
00209 zed3 - 3.556e-7*zed4;
00210 abund.solar[13] = (float)zsi;
00211
00212
00213 abund.solar[14] = abund.solar[13];
00214
00215
00216 zs = 1.12000;
00217 abund.solar[15] = (float)zs;
00218
00219
00220 zcl = 1.10000;
00221 abund.solar[16] = (float)zcl;
00222
00223
00224 zar = 1.091067724 + 2.51124e-6*zed3 - 0.0039589*sqrzed*zedlog +
00225 0.015686715*zedlog;
00226 abund.solar[17] = (float)zar;
00227
00228
00229 abund.solar[18] = abund.solar[13];
00230
00231
00232 zca = 1.077553875 - 0.00888806*zed + 0.001479866*zed2 - 6.5689e-5*
00233 zed3 + 1.16935e-6*zed4;
00234 abund.solar[19] = (float)zca;
00235
00236
00237 zfe = 0.223713045 + 0.001400746*zed + 0.000624734*zed2 - 3.5629e-5*
00238 zed3 + 8.13184e-7*zed4;
00239 abund.solar[25] = (float)zfe;
00240
00241
00242 abund.solar[20] = abund.solar[25];
00243
00244
00245 abund.solar[21] = abund.solar[25];
00246
00247
00248 abund.solar[22] = abund.solar[25];
00249
00250
00251 abund.solar[23] = abund.solar[25];
00252
00253
00254 abund.solar[24] = abund.solar[25];
00255
00256
00257 zco = zfe;
00258 abund.solar[26] = (float)zco;
00259
00260
00261 zni = zfe;
00262 abund.solar[27] = (float)zni;
00263
00264
00265 abund.solar[28] = abund.solar[25];
00266
00267
00268 abund.solar[29] = abund.solar[25];
00269
00270
00271 abund.solar[0] = 1.;
00272 abund.solar[1] = (float)(abund.solar[1]*abund.SolarSave[1]/
00273 zh);
00274
00275 for( i=2; i < LIMELM; i++ )
00276 {
00277 abund.solar[i] = (float)(abund.solar[i]*abund.SolarSave[i]*
00278 zed/zh);
00279 }
00280
00281 if( lgDebug )
00282 {
00283 fprintf( ioQQQ, "%10.2e", zed );
00284 for( i=0; i < LIMELM; i++ )
00285 {
00286 fprintf( ioQQQ, "%6.2f", log10(abund.solar[i]) );
00287 }
00288 fprintf( ioQQQ, "\n" );
00289
00290
00291 if( zed == zLimit )
00292
00293 {
00294 puts( "[Stop in abund_starburst]" );
00295 cdEXIT(EXIT_FAILURE);
00296 }
00297 }
00298
00299
00300 if( zed < zLimit )
00301 {
00302 zed = MIN2(zed*1.5,zLimit);
00303 }
00304 else
00305 {
00306 zed *= 1.5;
00307 }
00308 }
00309
00310
00311 if( optimize.lgVarOn )
00312 {
00313
00314 optimize.nvarxt[optimize.nparm] = 1;
00315 strcpy( optimize.chVarFmt[optimize.nparm], "ABUNDANCES STARBURST %f LOG" );
00316
00317 optimize.varang[optimize.nparm][0] = (float)log10(1e-3);
00318 optimize.varang[optimize.nparm][1] = (float)log10(zLimit);
00319
00320 optimize.nvfpnt[optimize.nparm] = input.nRead;
00321
00322 optimize.vparm[0][optimize.nparm] = (float)log10(zed);
00323 optimize.vincr[optimize.nparm] = 0.2f;
00324 ++optimize.nparm;
00325 }
00326
00327 DEBUG_EXIT( "abund_starburst()" );
00328 return;
00329 }
00330