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_sp[ipH_LIKE][nelem].numLevels_max;
00237 levels[ipHELIUM][nelem] = iso_sp[ipHE_LIKE][nelem].numLevels_max;
00238 iso_sp[ipH_LIKE][nelem].numLevels_max = 0;
00239 iso_sp[ipHE_LIKE][nelem].numLevels_max = 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 else if( p.nMatch(" ON ") )
00337 {
00338
00339 if( !lgElementSet )
00340 {
00341 fprintf( ioQQQ,
00342 " ParseElement did not find an element on the following line:\n" );
00343 p.PrintLine( ioQQQ );
00344 cdEXIT(EXIT_FAILURE);
00345 }
00346
00347 dense.lgElmtOn[nelem] = true;
00348
00349 if( levels[ipHYDROGEN][nelem] )
00350 {
00351 iso_sp[ipH_LIKE][nelem].numLevels_max = levels[ipHYDROGEN][nelem];
00352 iso_sp[ipHE_LIKE][nelem].numLevels_max = levels[ipHELIUM][nelem];
00353 }
00354 }
00355
00356 else if( p.nMatch("TABL") )
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 abund.lgAbunTabl[nelem] = true;
00368
00369
00370 abund.lgAbTaON = true;
00371
00372
00373
00374 if( p.nMatch("DEPT") )
00375 abund.lgAbTaDepth[nelem] = true;
00376 else
00377 abund.lgAbTaDepth[nelem] = false;
00378
00379
00380 if( nelem == ipHYDROGEN )
00381 {
00382 fprintf( ioQQQ, " cannot change abundance of hydrogen.\n" );
00383 fprintf( ioQQQ, " Sorry.\n" );
00384 cdEXIT(EXIT_FAILURE);
00385 }
00386
00387
00388 p.getline();
00389 abund.AbTabRad[0][nelem] = (realnum)p.FFmtRead();
00390 abund.AbTabFac[0][nelem] = (realnum)p.FFmtRead();
00391
00392 if( p.lgEOL() )
00393 {
00394 fprintf( ioQQQ, " no pairs entered - cannot interpolate\n" );
00395 cdEXIT(EXIT_FAILURE);
00396 }
00397
00398
00399 abund.nAbunTabl = 2;
00400
00401 lgEnd = false;
00402
00403
00404 while( !lgEnd && abund.nAbunTabl < LIMTABD )
00405 {
00406 p.getline();
00407 lgEnd = p.m_lgEOF;
00408 if( !lgEnd )
00409 {
00410
00411 if( p.strcmp("END") == 0 )
00412 lgEnd = true;
00413 }
00414
00415
00416 if( !lgEnd )
00417 {
00418 abund.AbTabRad[abund.nAbunTabl-1][nelem] =
00419 (realnum)p.FFmtRead();
00420
00421 abund.AbTabFac[abund.nAbunTabl-1][nelem] =
00422 (realnum)p.FFmtRead();
00423 abund.nAbunTabl += 1;
00424 }
00425 }
00426
00427 abund.nAbunTabl -= 1;
00428
00429
00430 for( long i=1; i < abund.nAbunTabl; i++ )
00431 {
00432
00433 if( abund.AbTabRad[i][nelem] <= abund.AbTabRad[i-1][nelem] )
00434 {
00435 fprintf( ioQQQ, "ParseElement: TABLE ELEMENT TABLE radii "
00436 "must be in increasing order\n" );
00437 cdEXIT(EXIT_FAILURE);
00438 }
00439 }
00440 }
00441
00442 else
00443 {
00444 fprintf( ioQQQ, "ParseElement: ELEMENT command - there must be "
00445 "a keyword on this line.\n" );
00446 fprintf( ioQQQ, " The keys I know about are TABLE, SCALE, _OFF, "
00447 "_ON_, IONIZATION, and ABUNDANCE.\n" );
00448 fprintf( ioQQQ, " Sorry.\n" );
00449 cdEXIT(EXIT_FAILURE);
00450 }
00451
00452
00453 if( optimize.lgVarOn )
00454 {
00455 if( p.nMatch("SCAL") )
00456 {
00457
00458 sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s SCALE %%f LOG",
00459 elementnames.chElementNameShort[nelem] );
00460
00461 }
00462 else if( p.nMatch("ABUN") )
00463 {
00464
00465 sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s ABUNDANCE %%f LOG",
00466 elementnames.chElementNameShort[nelem] );
00467 }
00468
00469 optimize.nvfpnt[optimize.nparm] = input.nRead;
00470 if( param >=0. )
00471 optimize.vparm[0][optimize.nparm] = (realnum)log10(param);
00472 else
00473 optimize.vparm[0][optimize.nparm] = (realnum)param;
00474 optimize.vincr[optimize.nparm] = 0.2f;
00475 optimize.nvarxt[optimize.nparm] = 1;
00476 ++optimize.nparm;
00477 }
00478 return;
00479 }