cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*abund_starburst generate abundance set from Fred Hamann's starburst evolution grid */ 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 /* this is largest stored metallicity */ 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 /* argument is metallicity */ 00057 zed = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00058 if( lgDebug ) 00059 { 00060 /* trace keyword did appear 00061 * this will not be used, but must trick the sanity test */ 00062 zHigh = zLimit; 00063 00064 /* if trace option set, print header now, and init ZED */ 00065 fprintf( ioQQQ, " Abundances relative to H, Z\n" ); 00066 00067 /* this is actual header line */ 00068 fprintf( ioQQQ, " Z "); 00069 00070 /* make line of chemical symbol names */ 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 /* no number and trace keyword did not appear */ 00082 fprintf( ioQQQ, " The metallicity must appear on this line.\n" ); 00083 cdEXIT(EXIT_FAILURE); 00084 } 00085 else 00086 { 00087 zHigh = zed; 00088 } 00089 00090 /* interpolate on Fred Hamann's set of starburst abundances 00091 * scan off keys _log or linear */ 00092 if( nMatch(" LOG",chCard) ) 00093 { 00094 zed = pow(10.,zed); 00095 zHigh = zed; 00096 } 00097 00098 else if( nMatch("LINE",chCard) ) 00099 { 00100 if( zed <= 0. ) 00101 { 00102 fprintf( ioQQQ, " Z .le.0 not allowed, Z=%10.2e\n", 00103 zed ); 00104 cdEXIT(EXIT_FAILURE); 00105 } 00106 } 00107 00108 else 00109 { 00110 /* log, linear not specified, neg so log */ 00111 if( zed <= 0. ) 00112 { 00113 zed = pow(10.,zed); 00114 zHigh = zed; 00115 } 00116 } 00117 00118 /* following if loop only if trace option is set 00119 * >>chng 97 jun 17, get rid of go to 00120 *999 continue 00121 * */ 00122 while( zed <= zHigh ) 00123 { 00124 if( zed < 1e-3 || zed > zLimit ) 00125 { 00126 fprintf( ioQQQ, " The metallicity must be between 1E-3 and%10.2e\n", 00127 zLimit ); 00128 cdEXIT(EXIT_FAILURE); 00129 } 00130 zed2 = zed*zed; 00131 zed3 = zed2*zed; 00132 zed4 = zed3*zed; 00133 zedlog = log(zed); 00134 sqrzed = sqrt(zed); 00135 /* The value of each element as a function of ZED=[Z] */ 00136 zh = 1.081646723 - 0.04600492*zed + 8.6569e-6*zed2 + 1.922e-5* 00137 zed3 - 2.0087e-7*zed4; 00138 00139 /* helium */ 00140 zhe = 0.864675891 + 0.044423807*zed + 7.10886e-5*zed2 - 5.3242e-5* 00141 zed3 + 5.70194e-7*zed4; 00142 abund.solar[1] = (realnum)zhe; 00143 00144 /* li, b, bo unchanged */ 00145 abund.solar[2] = 1.; 00146 abund.solar[3] = 1.; 00147 abund.solar[4] = 1.; 00148 00149 /* carbon */ 00150 zc = 0.347489799 + 0.016054107*zed - 0.00185855*zed2 + 5.43015e-5* 00151 zed3 - 5.3123e-7*zed4; 00152 abund.solar[5] = (realnum)zc; 00153 00154 /* nitrogen */ 00155 zn = -0.06549567 + 0.332073984*zed + 0.009146066*zed2 - 0.00054441* 00156 zed3 + 6.16363e-6*zed4; 00157 if( zn < 0.0 ) 00158 { 00159 zn = 0.05193*zed; 00160 } 00161 if( zed < 0.3 ) 00162 { 00163 zn = -0.00044731103 + 0.00026453554*zed + 0.52354843*zed2 - 00164 0.41156655*zed3 + 0.1290967*zed4; 00165 if( zn < 0.0 ) 00166 { 00167 zn = 0.000344828*zed; 00168 } 00169 } 00170 abund.solar[6] = (realnum)zn; 00171 00172 /* oxygen */ 00173 zo = 1.450212747 - 0.05379041*zed + 0.000498919*zed2 + 1.09646e-5* 00174 zed3 - 1.8147e-7*zed4; 00175 abund.solar[7] = (realnum)zo; 00176 00177 /* neon */ 00178 zne = 1.110139023 + 0.002551998*zed - 2.09516e-7*zed3 - 0.00798157* 00179 POW2(zedlog); 00180 abund.solar[9] = (realnum)zne; 00181 00182 /* fluorine, scale from neon */ 00183 abund.solar[8] = abund.solar[9]; 00184 00185 /* sodium */ 00186 zna = 1.072721387 - 0.02391599*POW2(zedlog) + .068644972* 00187 zedlog + 0.017030935/sqrzed; 00188 /* this one is slightly negative at very low Z */ 00189 zna = MAX2(1e-12,zna); 00190 abund.solar[10] = (realnum)zna; 00191 00192 /* magnesium */ 00193 zmg = 1.147209646 - 7.9491e-7*POW3(zed) - .00264458*POW2(zedlog) - 00194 0.00635552*zedlog; 00195 abund.solar[11] = (realnum)zmg; 00196 00197 /* aluminium */ 00198 zal = 1.068116358 - 0.00520227*sqrzed*zedlog - 0.01403851* 00199 POW2(zedlog) + 0.066186787*zedlog; 00200 /* this one is slightly negative at very low Z */ 00201 zal = MAX2(1e-12,zal); 00202 abund.solar[12] = (realnum)zal; 00203 00204 /* silicon */ 00205 zsi = 1.067372578 + 0.011818743*zed - 0.00107725*zed2 + 3.66056e-5* 00206 zed3 - 3.556e-7*zed4; 00207 abund.solar[13] = (realnum)zsi; 00208 00209 /* phosphorus scaled from silicon */ 00210 abund.solar[14] = abund.solar[13]; 00211 00212 /* sulphur */ 00213 zs = 1.12000; 00214 abund.solar[15] = (realnum)zs; 00215 00216 /* chlorine */ 00217 zcl = 1.10000; 00218 abund.solar[16] = (realnum)zcl; 00219 00220 /* argon */ 00221 zar = 1.091067724 + 2.51124e-6*zed3 - 0.0039589*sqrzed*zedlog + 00222 0.015686715*zedlog; 00223 abund.solar[17] = (realnum)zar; 00224 00225 /* potassium scaled from silicon */ 00226 abund.solar[18] = abund.solar[13]; 00227 00228 /* calcium */ 00229 zca = 1.077553875 - 0.00888806*zed + 0.001479866*zed2 - 6.5689e-5* 00230 zed3 + 1.16935e-6*zed4; 00231 abund.solar[19] = (realnum)zca; 00232 00233 /* iron */ 00234 zfe = 0.223713045 + 0.001400746*zed + 0.000624734*zed2 - 3.5629e-5* 00235 zed3 + 8.13184e-7*zed4; 00236 abund.solar[25] = (realnum)zfe; 00237 00238 /* scandium, scaled from iron */ 00239 abund.solar[20] = abund.solar[25]; 00240 00241 /* titanium, scaled from iron */ 00242 abund.solar[21] = abund.solar[25]; 00243 00244 /* vanadium scaled from iron */ 00245 abund.solar[22] = abund.solar[25]; 00246 00247 /* chromium scaled from iron */ 00248 abund.solar[23] = abund.solar[25]; 00249 00250 /* manganese scaled from iron */ 00251 abund.solar[24] = abund.solar[25]; 00252 00253 /* cobalt */ 00254 zco = zfe; 00255 abund.solar[26] = (realnum)zco; 00256 00257 /* nickel */ 00258 zni = zfe; 00259 abund.solar[27] = (realnum)zni; 00260 00261 /* copper scaled from iron */ 00262 abund.solar[28] = abund.solar[25]; 00263 00264 /* zinc scaled from iron */ 00265 abund.solar[29] = abund.solar[25]; 00266 00267 /* rescale to true abundances */ 00268 abund.solar[0] = 1.; 00269 abund.solar[1] = (realnum)(abund.solar[1]*abund.SolarSave[1]/ 00270 zh); 00271 00272 for( i=2; i < LIMELM; i++ ) 00273 { 00274 abund.solar[i] = (realnum)(abund.solar[i]*abund.SolarSave[i]* 00275 zed/zh); 00276 } 00277 00278 if( lgDebug ) 00279 { 00280 fprintf( ioQQQ, "%10.2e", zed ); 00281 for( i=0; i < LIMELM; i++ ) 00282 { 00283 fprintf( ioQQQ, "%6.2f", log10(abund.solar[i]) ); 00284 } 00285 fprintf( ioQQQ, "\n" ); 00286 00287 if( fp_equal( zed, zLimit ) ) 00288 { 00289 cdEXIT(EXIT_FAILURE); 00290 } 00291 } 00292 00293 /* this trick makes sure last one is upper limit */ 00294 if( zed < zLimit ) 00295 { 00296 zed = MIN2(zed*1.5,zLimit); 00297 } 00298 else 00299 { 00300 zed *= 1.5; 00301 } 00302 } 00303 00304 /* vary option */ 00305 if( optimize.lgVarOn ) 00306 { 00307 /* this is number of parameters to feed onto the input line */ 00308 optimize.nvarxt[optimize.nparm] = 1; 00309 strcpy( optimize.chVarFmt[optimize.nparm], "ABUNDANCES STARBURST %f LOG" ); 00310 /* following are upper and lower limits to metallicity range */ 00311 optimize.varang[optimize.nparm][0] = (realnum)log10(1e-3); 00312 optimize.varang[optimize.nparm][1] = (realnum)log10(zLimit); 00313 /* pointer to where to write */ 00314 optimize.nvfpnt[optimize.nparm] = input.nRead; 00315 /* log of metallicity will be argument */ 00316 optimize.vparm[0][optimize.nparm] = (realnum)log10(zed); 00317 optimize.vincr[optimize.nparm] = 0.2f; 00318 ++optimize.nparm; 00319 } 00320 return; 00321 }