00001
00002
00003
00004 #include "cddefines.h"
00005 #include "optimize.h"
00006 #include "abund.h"
00007 #include "dense.h"
00008 #include "input.h"
00009 #include "iso.h"
00010 #include "hmi.h"
00011 #include "mole.h"
00012 #include "elementnames.h"
00013 #include "parser.h"
00014
00015 void ParseElement( Parser &p)
00016 {
00017
00018
00019 static bool lgFirst = true;
00020 static long levels[NISO][LIMELM];
00021 bool lgEnd,
00022 lgHIT;
00023
00024 bool lgForceLog=false, lgForceLinear=false;
00025
00026 long int nelem,
00027 j,
00028 k;
00029 double param=0.0;
00030
00031 DEBUG_ENTRY( "ParseElement()" );
00032
00033
00034 if( lgFirst )
00035 {
00036 lgFirst = false;
00037 for( long i=0; i<NISO; ++i )
00038 {
00039 for( nelem=0; nelem<LIMELM; ++nelem )
00040 {
00041 levels[i][nelem] = 0;
00042 }
00043 }
00044 }
00045
00046 abund.lgAbnSolar = false;
00047
00048
00049 if( p.nMatch("READ") )
00050 {
00051 abund.npSolar = 0;
00052 for( long i=1; i < LIMELM; i++ )
00053 {
00054 p.getline();
00055
00056 if( p.m_lgEOF )
00057 {
00058 fprintf( ioQQQ, " Hit EOF while reading element list; use END to end list.\n" );
00059 cdEXIT(EXIT_FAILURE);
00060 }
00061
00062
00063 if( p.strcmp("END") == 0 )
00064 {
00065 return;
00066 }
00067
00068 j = 1;
00069 lgHIT = false;
00070 while( j < LIMELM && !lgHIT )
00071 {
00072 j += 1;
00073
00074 if( p.strcmp( elementnames.chElementNameShort[j-1] ) == 0 )
00075 {
00076 abund.npSolar += 1;
00077 abund.ipSolar[abund.npSolar-1] = j;
00078 lgHIT = true;
00079 }
00080 }
00081
00082 if( !lgHIT )
00083 {
00084 p.PrintLine( ioQQQ);
00085 fprintf( ioQQQ, " Sorry, but I did not recognize element name on this line.\n" );
00086 fprintf( ioQQQ, " Here is the list of names I recognize.\n" );
00087 fprintf( ioQQQ, " " );
00088
00089 for( k=2; k <= LIMELM; k++ )
00090 {
00091 fprintf( ioQQQ, "%4.4s\n", elementnames.chElementNameShort[k-1] );
00092 }
00093
00094 cdEXIT(EXIT_FAILURE);
00095 }
00096 }
00097
00098
00099 p.getline();
00100
00101 if( p.strcmp("END") == 0 )
00102 {
00103 return;
00104 }
00105
00106 else
00107 {
00108 fprintf( ioQQQ, " Too many elements were entered.\n" );
00109 fprintf( ioQQQ, " I only know about%3d elements.\n",
00110 LIMELM );
00111 fprintf( ioQQQ, " Sorry.\n" );
00112 cdEXIT(EXIT_FAILURE);
00113 }
00114 }
00115
00116
00117
00118 nelem = p.GetElem( );
00119
00120 bool lgElementSet = false;
00121 if( nelem >= ipHYDROGEN )
00122 {
00123 lgElementSet = true;
00124
00125 ASSERT( nelem>=0 && nelem < LIMELM );
00126 }
00127
00128
00129 if( p.nMatch(" LOG") )
00130 lgForceLog = true;
00131 else if( p.nMatch("LINE") )
00132 lgForceLinear = true;
00133
00134 if( p.nMatch("SCAL") )
00135 {
00136 param = p.FFmtRead();
00137
00138
00139 if( p.lgEOL() )
00140 {
00141 fprintf( ioQQQ, " There must be a number on this line.\n" );
00142 fprintf( ioQQQ, " Sorry.\n" );
00143 cdEXIT(EXIT_FAILURE);
00144 }
00145
00146 else
00147 {
00148 if( !lgElementSet )
00149 {
00150 fprintf( ioQQQ,
00151 " ParseElement did not find an element on the following line:\n" );
00152 p.PrintLine( ioQQQ );
00153 cdEXIT(EXIT_FAILURE);
00154 }
00155
00156 if( (lgForceLog || param <= 0.) && !lgForceLinear )
00157 {
00158
00159 param = pow(10.,param);
00160 }
00161 abund.ScaleElement[nelem] = (realnum)param;
00162 }
00163 }
00164
00165 else if( p.nMatch("ABUN") )
00166 {
00167 param = p.FFmtRead();
00168
00169
00170 if( p.lgEOL() )
00171 {
00172 fprintf( ioQQQ, " There must be a number on this line.\n" );
00173 fprintf( ioQQQ, " Sorry.\n" );
00174 cdEXIT(EXIT_FAILURE);
00175 }
00176
00177 else
00178 {
00179 if( !lgElementSet )
00180 {
00181 fprintf( ioQQQ,
00182 " ParseElement did not find an element on the following line:\n" );
00183 p.PrintLine( ioQQQ );
00184 cdEXIT(EXIT_FAILURE);
00185 }
00186
00187 if( !lgForceLinear )
00188 {
00189
00190 param = pow(10.,param);
00191 }
00192 abund.solar[nelem] = (realnum)param;
00193
00194 if( abund.solar[nelem] > 1. )
00195 {
00196 fprintf( ioQQQ,
00197 " Please check the abundance of this element. It seems high to me.\n" );
00198 }
00199 }
00200 }
00201
00202 else if( p.nMatch(" OFF") )
00203 {
00204 param = p.FFmtRead();
00205
00206
00207
00208
00209 if( p.nMatch( "LIMI") )
00210 {
00211 if( p.lgEOL() )
00212 {
00213 fprintf( ioQQQ, " There must be an abundances on the ELEMENT OFF LIMIT command.\n" );
00214 fprintf( ioQQQ, " Sorry.\n" );
00215 cdEXIT(EXIT_FAILURE);
00216 }
00217 dense.AbundanceLimit = (realnum)pow(10., param);
00218 }
00219 else
00220 {
00221 if( !lgElementSet )
00222 {
00223 fprintf( ioQQQ,
00224 " ParseElement did not find an element on the following line:\n" );
00225 p.PrintLine ( ioQQQ );
00226 cdEXIT(EXIT_FAILURE);
00227 }
00228
00229 dense.lgElmtOn[nelem] = false;
00230
00231 dense.IonHigh[nelem] = -1;
00232 dense.lgSetIoniz[nelem] = false;
00233
00234 if( nelem > ipHELIUM )
00235 {
00236 levels[ipHYDROGEN][nelem] = iso.numLevels_max[ipH_LIKE][nelem];
00237 levels[ipHELIUM][nelem] = iso.numLevels_max[ipHE_LIKE][nelem];
00238 iso.numLevels_max[ipH_LIKE][nelem] = 0;
00239 iso.numLevels_max[ipHE_LIKE][nelem] = 0;
00240 }
00241
00242 if( nelem == ipHYDROGEN )
00243 {
00244 fprintf( ioQQQ, " It is not possible to turn hydrogen off.\n" );
00245 fprintf( ioQQQ, " Sorry.\n" );
00246 cdEXIT(EXIT_FAILURE);
00247 }
00248 }
00249 }
00250
00251
00252 else if( p.nMatch("IONI") )
00253 {
00254 bool lgLogSet = false;
00255 int ion;
00256 int ihi , low;
00257
00258
00259 if( !lgElementSet )
00260 {
00261 fprintf( ioQQQ,
00262 " ParseElement did not find an element on the following line:\n" );
00263 p.PrintLine( ioQQQ );
00264 cdEXIT(EXIT_FAILURE);
00265 }
00266 if( !dense.lgElmtOn[nelem] )
00267 {
00268 fprintf(ioQQQ,"Sorry, you cannot set the ionization of %s since it has been turned off.\n" ,
00269 elementnames.chElementName[nelem] );
00270 cdEXIT(EXIT_FAILURE);
00271 }
00272
00273
00274 dense.lgSetIoniz[nelem] = true;
00275 ion = 0;
00276 while( ion<nelem+2 )
00277 {
00278
00279
00280
00281 dense.SetIoniz[nelem][ion] = (realnum)p.FFmtRead();
00282
00283 if( p.lgEOL() )
00284 break;
00285
00286
00287 if( dense.SetIoniz[nelem][ion] < 0. )
00288 lgLogSet = true;
00289 ++ion;
00290 }
00291 param = dense.SetIoniz[nelem][0];
00292
00293
00294 for( long i=ion; i<nelem+2; ++i )
00295 dense.SetIoniz[nelem][i] = 0.;
00296
00297
00298 if( lgLogSet )
00299 {
00300 for( long i=0; i<ion; ++i )
00301 dense.SetIoniz[nelem][i] = (realnum)pow((realnum)10.f, dense.SetIoniz[nelem][i]);
00302 }
00303
00304
00305 low = 0;
00306 while( dense.SetIoniz[nelem][low]==0 && low < nelem+1 )
00307 ++low;
00308 ihi = nelem+1;
00309 while( dense.SetIoniz[nelem][ihi]==0 && ihi > low )
00310 --ihi;
00311
00312 if( ihi==low && dense.SetIoniz[nelem][ihi]==0 )
00313 {
00314 fprintf(ioQQQ," element ionization command has all zero ionization fractions. This is not possible.\n Sorry\n");
00315 cdEXIT(EXIT_FAILURE);
00316 }
00317 realnum AbundSum=0.;
00318 for(long i=low; i<=ihi; ++i )
00319 {
00320 if( dense.SetIoniz[nelem][i]==0 )
00321 {
00322 fprintf(ioQQQ," element abundance command has zero abundance between positive values. This is not possible.\n Sorry\n");
00323 cdEXIT(EXIT_FAILURE);
00324 }
00325 AbundSum += dense.SetIoniz[nelem][i];
00326 }
00327
00328
00329 for(long i=low; i<=ihi; ++i )
00330 dense.SetIoniz[nelem][i] /= AbundSum;
00331
00332
00333
00334
00335
00336
00337 if( nelem==ipHYDROGEN && !hmi.lgNoH2Mole )
00338 {
00339 fprintf(ioQQQ," Hydrogen ionization has been set, H chemistry is disabled.\n");
00340
00341 hmi.lgNoH2Mole = true;
00342 }
00343
00344
00345
00346 if( mole.lgElem_in_chemistry[nelem] && !co.lgNoCOMole )
00347 {
00348 fprintf(ioQQQ," The ionization of an element in the CO chemistry network has been set, CO chemistry is disabled.\n");
00349
00350 co.lgNoCOMole = true;
00351 }
00352 }
00353
00354
00355
00356 else if( p.nMatch(" ON ") )
00357 {
00358
00359 if( !lgElementSet )
00360 {
00361 fprintf( ioQQQ,
00362 " ParseElement did not find an element on the following line:\n" );
00363 p.PrintLine( ioQQQ );
00364 cdEXIT(EXIT_FAILURE);
00365 }
00366
00367 dense.lgElmtOn[nelem] = true;
00368
00369 if( levels[ipHYDROGEN][nelem] )
00370 {
00371 iso.numLevels_max[ipH_LIKE][nelem] = levels[ipHYDROGEN][nelem];
00372 iso.numLevels_max[ipHE_LIKE][nelem] = levels[ipHELIUM][nelem];
00373 }
00374 }
00375
00376 else if( p.nMatch("TABL") )
00377 {
00378
00379 if( !lgElementSet )
00380 {
00381 fprintf( ioQQQ,
00382 " ParseElement did not find an element on the following line:\n" );
00383 p.PrintLine( ioQQQ );
00384 cdEXIT(EXIT_FAILURE);
00385 }
00386
00387 abund.lgAbunTabl[nelem] = true;
00388
00389
00390 abund.lgAbTaON = true;
00391
00392
00393
00394 if( p.nMatch("DEPT") )
00395 abund.lgAbTaDepth[nelem] = true;
00396 else
00397 abund.lgAbTaDepth[nelem] = false;
00398
00399
00400 if( nelem == ipHYDROGEN )
00401 {
00402 fprintf( ioQQQ, " cannot change abundance of hydrogen.\n" );
00403 fprintf( ioQQQ, " Sorry.\n" );
00404 cdEXIT(EXIT_FAILURE);
00405 }
00406
00407
00408 p.getline();
00409 abund.AbTabRad[0][nelem] = (realnum)p.FFmtRead();
00410 abund.AbTabFac[0][nelem] = (realnum)p.FFmtRead();
00411
00412 if( p.lgEOL() )
00413 {
00414 fprintf( ioQQQ, " no pairs entered - cannot interpolate\n" );
00415 cdEXIT(EXIT_FAILURE);
00416 }
00417
00418
00419 abund.nAbunTabl = 2;
00420
00421 lgEnd = false;
00422
00423
00424 while( !lgEnd && abund.nAbunTabl < LIMTABD )
00425 {
00426 p.getline();
00427 lgEnd = p.m_lgEOF;
00428 if( !lgEnd )
00429 {
00430
00431 if( p.strcmp("END") == 0 )
00432 lgEnd = true;
00433 }
00434
00435
00436 if( !lgEnd )
00437 {
00438 abund.AbTabRad[abund.nAbunTabl-1][nelem] =
00439 (realnum)p.FFmtRead();
00440
00441 abund.AbTabFac[abund.nAbunTabl-1][nelem] =
00442 (realnum)p.FFmtRead();
00443 abund.nAbunTabl += 1;
00444 }
00445 }
00446
00447 abund.nAbunTabl -= 1;
00448
00449
00450 for( long i=1; i < abund.nAbunTabl; i++ )
00451 {
00452
00453 if( abund.AbTabRad[i][nelem] <= abund.AbTabRad[i-1][nelem] )
00454 {
00455 fprintf( ioQQQ, "ParseElement: TABLE ELEMENT TABLE radii "
00456 "must be in increasing order\n" );
00457 cdEXIT(EXIT_FAILURE);
00458 }
00459 }
00460 }
00461
00462 else
00463 {
00464 fprintf( ioQQQ, "ParseElement: ELEMENT command - there must be "
00465 "a keyword on this line.\n" );
00466 fprintf( ioQQQ, " The keys I know about are TABLE, SCALE, _OFF, "
00467 "_ON_, IONIZATION, and ABUNDANCE.\n" );
00468 fprintf( ioQQQ, " Sorry.\n" );
00469 cdEXIT(EXIT_FAILURE);
00470 }
00471
00472
00473 if( optimize.lgVarOn )
00474 {
00475 if( p.nMatch("SCAL") )
00476 {
00477
00478 sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s SCALE %%f LOG",
00479 elementnames.chElementNameShort[nelem] );
00480
00481 }
00482 else if( p.nMatch("ABUN") )
00483 {
00484
00485 sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s ABUNDANCE %%f LOG",
00486 elementnames.chElementNameShort[nelem] );
00487 }
00488
00489 optimize.nvfpnt[optimize.nparm] = input.nRead;
00490 if( param >=0. )
00491 optimize.vparm[0][optimize.nparm] = (realnum)log10(param);
00492 else
00493 optimize.vparm[0][optimize.nparm] = (realnum)param;
00494 optimize.vincr[optimize.nparm] = 0.2f;
00495 optimize.nvarxt[optimize.nparm] = 1;
00496 ++optimize.nparm;
00497 }
00498 return;
00499 }