cloudy: /home66/gary/public_html/cloudy/c08_branch/source/abund

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 }
/home66/gary/public_html/cloudy/c08_branch/source/abund_starburst.cpp
