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