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 /*ParseElement parse options on element command */ 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 "parse.h" 00014 00015 void ParseElement(char *chCard ) 00016 { 00017 /* this will be used to remember how many levels are active in any element we turn off, 00018 * so that we can retain state if turned back on */ 00019 static bool lgFirst = true; 00020 static long levels[NISO][LIMELM]; 00021 char chCap[INPUT_LINE_LENGTH]; 00022 bool lgEOL, 00023 lgEnd, 00024 lgHIT; 00025 00026 bool lgForceLog=false, lgForceLinear=false; 00027 00028 long int i, 00029 nelem, 00030 j, 00031 k; 00032 double param; 00033 00034 DEBUG_ENTRY( "ParseElement()" ); 00035 00036 /* zero out this array if first call */ 00037 if( lgFirst ) 00038 { 00039 lgFirst = false; 00040 for( i=0; i<NISO; ++i ) 00041 { 00042 for( nelem=0; nelem<LIMELM; ++nelem ) 00043 { 00044 levels[i][nelem] = 0; 00045 } 00046 } 00047 } 00048 /* say that abundances have been changed */ 00049 abund.lgAbnSolar = false; 00050 00051 /* read in list of elements for abundances command */ 00052 if( nMatch("READ",chCard) ) 00053 { 00054 abund.npSolar = 0; 00055 for( i=1; i < LIMELM; i++ ) 00056 { 00057 input_readarray(chCard,&lgEnd); 00058 00059 if( lgEnd ) 00060 { 00061 fprintf( ioQQQ, " Hit EOF while reading element list; use END to end list.\n" ); 00062 cdEXIT(EXIT_FAILURE); 00063 } 00064 00065 strcpy( chCap, chCard ); 00066 caps(chCap); 00067 00068 /* this would be a line starting with END to say end on list */ 00069 if( strncmp(chCap,"END",3) == 0 ) 00070 { 00071 return; 00072 } 00073 00074 j = 1; 00075 lgHIT = false; 00076 while( j < LIMELM && !lgHIT ) 00077 { 00078 j += 1; 00079 00080 if( strncmp( chCap,elementnames.chElementNameShort[j-1] , 4) == 0 ) 00081 { 00082 abund.npSolar += 1; 00083 abund.ipSolar[abund.npSolar-1] = j; 00084 lgHIT = true; 00085 } 00086 } 00087 00088 if( !lgHIT ) 00089 { 00090 fprintf( ioQQQ, "%80.80s\n", chCard ); 00091 fprintf( ioQQQ, " Sorry, but I did not recognize element name on this line.\n" ); 00092 fprintf( ioQQQ, " Here is the list of names I recognize.\n" ); 00093 fprintf( ioQQQ, " " ); 00094 00095 for( k=2; k <= LIMELM; k++ ) 00096 { 00097 fprintf( ioQQQ, "%4.4s\n", elementnames.chElementNameShort[k-1] ); 00098 } 00099 00100 cdEXIT(EXIT_FAILURE); 00101 } 00102 } 00103 00104 /* fell through, make sure one more line, the end line */ 00105 input_readarray(chCard,&lgEnd); 00106 strcpy( chCap, chCard ); 00107 caps(chCap); 00108 00109 if( strncmp(chCap,"END",3) == 0 ) 00110 { 00111 return; 00112 } 00113 00114 else 00115 { 00116 fprintf( ioQQQ, " Too many elements were entered.\n" ); 00117 fprintf( ioQQQ, " I only know about%3d elements.\n", 00118 LIMELM ); 00119 fprintf( ioQQQ, " Sorry.\n" ); 00120 cdEXIT(EXIT_FAILURE); 00121 } 00122 } 00123 00124 /* find which element - will be used in remainder of routine to 00125 * adjust aspects of this element */ 00126 nelem = GetElem(chCard ); 00127 00128 bool lgElementSet = false; 00129 if( nelem >= ipHYDROGEN ) 00130 { 00131 lgElementSet = true; 00132 /* nelem is now the atomic number of the element, and must be correct */ 00133 ASSERT( nelem>=0 && nelem < LIMELM ); 00134 } 00135 00136 /* look for log or linear keywords */ 00137 if( nMatch(" LOG",chCard) ) 00138 lgForceLog = true; 00139 else if( nMatch("LINE",chCard) ) 00140 lgForceLinear = true; 00141 00142 i = 4; 00143 param = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00144 00145 if( nMatch("SCAL",chCard) ) 00146 { 00147 /* enter abundance as scale factor, relative to what is in now */ 00148 if( lgEOL ) 00149 { 00150 fprintf( ioQQQ, " There must be a number on this line.\n" ); 00151 fprintf( ioQQQ, " Sorry.\n" ); 00152 cdEXIT(EXIT_FAILURE); 00153 } 00154 00155 else 00156 { 00157 00158 if( !lgElementSet ) 00159 { 00160 fprintf( ioQQQ, 00161 " ParseElement did not find an element on the following line:\n" ); 00162 fprintf( ioQQQ, 00163 "%80.80s\n", chCard ); 00164 cdEXIT(EXIT_FAILURE); 00165 } 00166 /* interpret as log unless forced linear */ 00167 if( (lgForceLog || param <= 0.) && !lgForceLinear ) 00168 { 00169 /* option for log of scale factor */ 00170 param = pow(10.,param); 00171 } 00172 abund.ScaleElement[nelem] = (realnum)param; 00173 } 00174 } 00175 00176 else if( nMatch("ABUN",chCard) ) 00177 { 00178 /* log of actual abundance */ 00179 if( lgEOL ) 00180 { 00181 fprintf( ioQQQ, " There must be a number on this line.\n" ); 00182 fprintf( ioQQQ, " Sorry.\n" ); 00183 cdEXIT(EXIT_FAILURE); 00184 } 00185 00186 else 00187 { 00188 00189 if( !lgElementSet ) 00190 { 00191 fprintf( ioQQQ, 00192 " ParseElement did not find an element on the following line:\n" ); 00193 fprintf( ioQQQ, 00194 "%80.80s\n", chCard ); 00195 cdEXIT(EXIT_FAILURE); 00196 } 00197 if( lgForceLinear ) 00198 { 00199 abund.solar[nelem] = (realnum)param; 00200 } 00201 else 00202 { 00203 abund.solar[nelem] = (realnum)pow(10.,param); 00204 } 00205 00206 if( abund.solar[nelem]> 1. ) 00207 { 00208 fprintf( ioQQQ, 00209 " Please check the abundance of this element. It seems high to me.\n" ); 00210 } 00211 } 00212 } 00213 00214 else if( nMatch(" OFF",chCard) ) 00215 { 00216 /* keyword LIMIT sets lower limit to abundance of element 00217 * that will be included */ 00218 /* option to turn off this element, set abundance to zero */ 00219 if( nMatch( "LIMI",chCard) ) 00220 { 00221 if( lgEOL ) 00222 { 00223 fprintf( ioQQQ, " There must be an abundances on the ELEMENT OFF LIMIT command.\n" ); 00224 fprintf( ioQQQ, " Sorry.\n" ); 00225 cdEXIT(EXIT_FAILURE); 00226 } 00227 dense.AbundanceLimit = (realnum)pow(10., param); 00228 } 00229 else 00230 { 00231 if( !lgElementSet ) 00232 { 00233 fprintf( ioQQQ, 00234 " ParseElement did not find an element on the following line:\n" ); 00235 fprintf( ioQQQ, 00236 "%80.80s\n", chCard ); 00237 cdEXIT(EXIT_FAILURE); 00238 } 00239 /* specify element name */ 00240 dense.lgElmtOn[nelem] = false; 00241 /* another flag that element is off */ 00242 dense.IonHigh[nelem] = -1; 00243 dense.lgSetIoniz[nelem] = false; 00244 /* set no levels for all elements that are turned off, except for helium itself, which always exists */ 00245 if( nelem > ipHELIUM ) 00246 { 00247 levels[ipHYDROGEN][nelem] = iso.numLevels_max[ipH_LIKE][nelem]; 00248 levels[ipHELIUM][nelem] = iso.numLevels_max[ipHE_LIKE][nelem]; 00249 iso.numLevels_max[ipH_LIKE][nelem] = 0; 00250 iso.numLevels_max[ipHE_LIKE][nelem] = 0; 00251 } 00252 00253 if( nelem == ipHYDROGEN ) 00254 { 00255 fprintf( ioQQQ, " It is not possible to turn hydrogen off.\n" ); 00256 fprintf( ioQQQ, " Sorry.\n" ); 00257 cdEXIT(EXIT_FAILURE); 00258 } 00259 } 00260 } 00261 00262 /* specify an ionization distribution */ 00263 else if( nMatch("IONI",chCard) ) 00264 { 00265 bool lgLogSet = false; 00266 int ion; 00267 int ihi , low; 00268 00269 00270 if( !lgElementSet ) 00271 { 00272 fprintf( ioQQQ, 00273 " ParseElement did not find an element on the following line:\n" ); 00274 fprintf( ioQQQ, 00275 "%80.80s\n", chCard ); 00276 cdEXIT(EXIT_FAILURE); 00277 } 00278 if( !dense.lgElmtOn[nelem] ) 00279 { 00280 fprintf(ioQQQ,"Sorry, you cannot set the ionization of %s since it has been turned off.\n" , 00281 elementnames.chElementName[nelem] ); 00282 cdEXIT(EXIT_FAILURE); 00283 } 00284 00285 i = 4; 00286 /* option to specify ionization distribution for this element */ 00287 dense.lgSetIoniz[nelem] = true; 00288 ion = 0; 00289 while( ion<nelem+2 ) 00290 { 00291 /* the ionization fractions that are set when above set true, 00292 * gas phase abundance is this times total abundance 00293 * Ionization fraction for [nelem][ion] */ 00294 dense.SetIoniz[nelem][ion] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00295 00296 if( lgEOL ) 00297 break; 00298 00299 /* all are log if any are negative */ 00300 if( dense.SetIoniz[nelem][ion] < 0. ) 00301 lgLogSet = true; 00302 ++ion; 00303 } 00304 00305 /* zero out ones not specified */ 00306 for( i=ion; i<nelem+2; ++i ) 00307 dense.SetIoniz[nelem][i] = 0.; 00308 00309 /* convert rest to linear if any were negative */ 00310 if( lgLogSet ) 00311 { 00312 for( i=0; i<ion; ++i ) 00313 dense.SetIoniz[nelem][i] = (realnum)pow((realnum)10.f, dense.SetIoniz[nelem][i]); 00314 } 00315 00316 /* now check that no zero abundances exist between lowest and highest non-zero values */ 00317 low = 0; 00318 while( dense.SetIoniz[nelem][low]==0 && low < nelem+1 ) 00319 ++low; 00320 ihi = nelem+1; 00321 while( dense.SetIoniz[nelem][ihi]==0 && ihi > low ) 00322 --ihi; 00323 00324 if( ihi==low && dense.SetIoniz[nelem][ihi]==0 ) 00325 { 00326 fprintf(ioQQQ," element ionization command has all zero ionization fractions. This is not possible.\n Sorry\n"); 00327 cdEXIT(EXIT_FAILURE); 00328 } 00329 for(i=low; i<=ihi; ++i ) 00330 { 00331 if( dense.SetIoniz[nelem][i]==0 ) 00332 { 00333 fprintf(ioQQQ," element abundance command has zero abundance between positive values. This is not possible.\n Sorry\n"); 00334 cdEXIT(EXIT_FAILURE); 00335 } 00336 } 00337 /* >>chng 05 mar 13, add this 00338 * finally turn off any part of the chemical networks that might use this element, 00339 * since chemistry must come into equilibrium with the ionization, and vice vese - 00340 * ion distribution cannot come into equilibrium with chemistry when ions are set, 00341 * so much turn off chemistry */ 00342 if( nelem==ipHYDROGEN && !hmi.lgNoH2Mole ) 00343 { 00344 fprintf(ioQQQ," Hydrogen ionization has been set, H chemistry is disabled.\n"); 00345 /* turn off only H2 */ 00346 hmi.lgNoH2Mole = true; 00347 } 00348 /* loop over all atoms in co chem net - these occur as the last NUM_ELEMENTS 00349 * elements in the full set of mole.num_heavy_molec */ 00350 /* >>chng 06 mar 19, do not print if network already turned off */ 00351 if( mole.lgElem_in_chemistry[nelem] && !co.lgNoCOMole ) 00352 { 00353 fprintf(ioQQQ," The ionization of an element in the CO chemistry network has been set, CO chemistry is disabled.\n"); 00354 /* turn off CO network */ 00355 co.lgNoCOMole = true; 00356 } 00357 } 00358 00359 /* option to turn element on - some ini files turn off elements, 00360 * this can turn them back on */ 00361 else if( nMatch(" ON ",chCard) ) 00362 { 00363 00364 if( !lgElementSet ) 00365 { 00366 fprintf( ioQQQ, 00367 " ParseElement did not find an element on the following line:\n" ); 00368 fprintf( ioQQQ, 00369 "%80.80s\n", chCard ); 00370 cdEXIT(EXIT_FAILURE); 00371 } 00372 /* option to turn off this element, set abundance to zero */ 00373 dense.lgElmtOn[nelem] = true; 00374 /* reset levels to default if they were ever turned off with element off command */ 00375 if( levels[ipHYDROGEN][nelem] ) 00376 { 00377 iso.numLevels_max[ipH_LIKE][nelem] = levels[ipHYDROGEN][nelem]; 00378 iso.numLevels_max[ipHE_LIKE][nelem] = levels[ipHELIUM][nelem]; 00379 } 00380 } 00381 00382 else if( nMatch("TABL",chCard) ) 00383 { 00384 00385 if( !lgElementSet ) 00386 { 00387 fprintf( ioQQQ, 00388 " ParseElement did not find an element on the following line:\n" ); 00389 fprintf( ioQQQ, 00390 "%80.80s\n", chCard ); 00391 cdEXIT(EXIT_FAILURE); 00392 } 00393 /* read in table of position-dependent abundances for a particular element. */ 00394 abund.lgAbunTabl[nelem] = true; 00395 00396 /* general flag saying this option turned on */ 00397 abund.lgAbTaON = true; 00398 00399 /* does the table give depth or radius? keyword DEPTH 00400 * says it is depth, default is radius */ 00401 if( nMatch("DEPT",chCard) ) 00402 abund.lgAbTaDepth[nelem] = true; 00403 else 00404 abund.lgAbTaDepth[nelem] = false; 00405 00406 /* make sure not trying to change hydrogen */ 00407 if( nelem == ipHYDROGEN ) 00408 { 00409 fprintf( ioQQQ, " cannot change abundance of hydrogen.\n" ); 00410 fprintf( ioQQQ, " Sorry.\n" ); 00411 cdEXIT(EXIT_FAILURE); 00412 } 00413 00414 /* read pair giving radius and abundance */ 00415 input_readarray(chCard,&lgEnd); 00416 i = 1; 00417 abund.AbTabRad[0][nelem] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00418 abund.AbTabFac[0][nelem] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00419 00420 if( lgEOL ) 00421 { 00422 fprintf( ioQQQ, " no pairs entered - cannot interpolate\n" ); 00423 cdEXIT(EXIT_FAILURE); 00424 } 00425 00426 /* number of points in the table */ 00427 abund.nAbunTabl = 2; 00428 00429 lgEnd = false; 00430 /* LIMTAB is limit to number of points we can store and is 00431 * set to 500 in abundances */ 00432 while( !lgEnd && abund.nAbunTabl < LIMTABD ) 00433 { 00434 input_readarray(chCard,&lgEnd); 00435 if( !lgEnd ) 00436 { 00437 /* convert first 4 char to caps, into chCap */ 00438 cap4(chCap , chCard); 00439 if( strncmp(chCap,"END",3) == 0 ) 00440 lgEnd = true; 00441 } 00442 00443 /* lgEnd may have been set within above if, if end line encountered*/ 00444 if( !lgEnd ) 00445 { 00446 i = 1; 00447 abund.AbTabRad[abund.nAbunTabl-1][nelem] = 00448 (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL); 00449 00450 abund.AbTabFac[abund.nAbunTabl-1][nelem] = 00451 (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL); 00452 abund.nAbunTabl += 1; 00453 } 00454 } 00455 00456 abund.nAbunTabl -= 1; 00457 00458 /* now check that radii are in increasing order */ 00459 for( i=1; i < abund.nAbunTabl; i++ ) 00460 { 00461 /* the radius values are assumed to be strictly increasing */ 00462 if( abund.AbTabRad[i][nelem] <= abund.AbTabRad[i-1][nelem] ) 00463 { 00464 fprintf( ioQQQ, "ParseElement: TABLE ELEMENT TABLE radii " 00465 "must be in increasing order\n" ); 00466 cdEXIT(EXIT_FAILURE); 00467 } 00468 } 00469 } 00470 00471 else 00472 { 00473 fprintf( ioQQQ, "ParseElement: ELEMENT command - there must be " 00474 "a keyword on this line.\n" ); 00475 fprintf( ioQQQ, " The keys I know about are TABLE, SCALE, _OFF, " 00476 "_ON_, IONIZATION, and ABUNDANCE.\n" ); 00477 fprintf( ioQQQ, " Sorry.\n" ); 00478 cdEXIT(EXIT_FAILURE); 00479 } 00480 00481 /* vary option */ 00482 if( optimize.lgVarOn ) 00483 { 00484 optimize.nvarxt[optimize.nparm] = 1; 00485 /* pointer to where to write */ 00486 optimize.nvfpnt[optimize.nparm] = input.nRead; 00487 00488 if( nMatch("SCAL",chCard) ) 00489 { 00490 /* vary scale factor */ 00491 sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s SCALE %%f LOG", 00492 elementnames.chElementNameShort[nelem] ); 00493 00494 /* param is linear scale factor */ 00495 optimize.vparm[0][optimize.nparm] = (realnum)log10(param); 00496 optimize.vincr[optimize.nparm] = 0.2f; 00497 } 00498 00499 else if( nMatch("ABUN",chCard) ) 00500 { 00501 /* vary absolute abundance */ 00502 sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s ABUND %%f LOG ", 00503 elementnames.chElementNameShort[nelem] ); 00504 00505 /* param is log of abundance by number relative to hydrogen */ 00506 optimize.vparm[0][optimize.nparm] = (realnum)param; 00507 optimize.vincr[optimize.nparm] = 0.2f; 00508 } 00509 ++optimize.nparm; 00510 } 00511 return; 00512 }