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 /*ParseAtomISO parse information from the atom XX-like command line */ 00004 #include "cddefines.h" 00005 #include "elementnames.h" 00006 #include "helike_recom.h" 00007 #include "hydrogenic.h" 00008 #include "iso.h" 00009 #include "parse.h" 00010 #include "phycon.h" 00011 #include "rfield.h" 00012 #include "taulines.h" 00013 #include "thirdparty.h" 00014 00015 /*ParseAtomISO parse parameters off the XX-like command */ 00016 void ParseAtomISO(long ipISO, char *chCard ) 00017 { 00018 bool lgEOL; 00019 long int i, 00020 numLevels, 00021 nelem; 00022 00023 DEBUG_ENTRY( "ParseAtomISO()" ); 00024 00025 /* look for the name of an element - if we don't find one do the entire 00026 * iso sequence - returns negative number if element not found */ 00027 nelem = GetElem( chCard ); 00028 00029 /* H-like Helium is not possible */ 00030 if( ipISO==ipHE_LIKE && nelem==ipHYDROGEN ) 00031 { 00032 fprintf(ioQQQ," Sorry, H-like He does not exist.\n"); 00033 cdEXIT(EXIT_FAILURE); 00034 } 00035 00036 /* collisions - don't pick up the levels collapsed command */ 00037 if( nMatch("COLL",chCard) && !nMatch("LEVE" , chCard ) ) 00038 { 00039 /* option to turn collisions off, all are 00040 * set to 1 in zero. command can accept only one option at a time */ 00041 if( nMatch("EXCI",chCard) ) 00042 { 00043 /* turn off collisional excitation */ 00044 iso.lgColl_excite[ipISO] = false; 00045 phycon.lgPhysOK = false; 00046 } 00047 else if( nMatch("IONI",chCard) ) 00048 { 00049 /* turn off collisional ionization */ 00050 iso.lgColl_ionize[ipISO] = false; 00051 phycon.lgPhysOK = false; 00052 } 00053 00054 else if( nMatch("2S2P",chCard) || ( nMatch("2P2S",chCard) && ipISO == ipH_LIKE ) ) 00055 { 00056 /* >>chng 02 feb 07, change from 2s2p to l-mixing */ 00057 /* this is the atom h-like collision l-mixing command */ 00058 fprintf(ioQQQ,"This command changed to ATOM H-LIKE COLLISIONS L-MIXING\n"); 00059 fprintf(ioQQQ,"I will parse it for now, but may not in the future.\n"); 00060 /* turn off 2s - 2p collisions */ 00061 iso.lgColl_l_mixing[ipISO] = false; 00062 phycon.lgPhysOK = false; 00063 } 00064 00065 else if( nMatch("L-MI",chCard) ) 00066 { 00067 if( ipISO == ipH_LIKE ) 00068 { 00069 /* >>chng 02 feb 07, change from 2s2p to l-mixing */ 00070 /* this is the atom h-like collision l-mixing command */ 00071 iso.lgColl_l_mixing[ipISO] = false; 00072 phycon.lgPhysOK = false; 00073 } 00074 else if( nMatch("THER",chCard) ) 00075 { 00076 /* use l-mix from 00077 *>>refer l-mix all Vrinceanu, D. & Flannery, M. R. 2001, PhysRevA 63, 032701 */ 00078 if( nMatch("NO T",chCard) ) 00079 { 00080 /* This is the "NO Thermal average" command. It 00081 * causes collisions strengths to be evaluated at kT rather than 00082 * integrated over a Maxwellian peaked at kT. */ 00083 iso.lgCS_therm_ave[ipISO] = false; 00084 } 00085 else 00086 { 00087 iso.lgCS_therm_ave[ipISO] = true; 00088 } 00089 } 00090 else if( nMatch("PENG",chCard) ) 00091 { 00092 iso.lgCS_Vrinceanu[ipISO] = false; 00093 } 00094 else if( nMatch(" OFF" , chCard ) ) 00095 { 00096 /* this is the atom xx-like collision l-mixing command */ 00097 /* turn off same-n collisions */ 00098 iso.lgColl_l_mixing[ipISO] = false; 00099 phycon.lgPhysOK = false; 00100 iso.lgCS_Vrinceanu[ipISO] = false; 00101 iso.lgCS_therm_ave[ipISO] = false; 00102 } 00103 else 00104 { 00105 fprintf( ioQQQ, " needs parameter\n" ); 00106 cdEXIT(EXIT_FAILURE); 00107 } 00108 } 00109 else if( nMatch(" OFF" , chCard ) ) 00110 { 00111 /* turn everything off, since no keyword given */ 00112 iso.lgColl_excite[ipISO] = false; 00113 iso.lgColl_ionize[ipISO] = false; 00114 iso.lgColl_l_mixing[ipISO] = false; 00115 phycon.lgPhysOK = false; 00116 } 00117 else 00118 { 00119 fprintf( ioQQQ, " needs parameter\n" ); 00120 cdEXIT(EXIT_FAILURE); 00121 } 00122 } 00123 00124 else if( nMatch("DAMP",chCard) ) 00125 { 00126 if( ipISO == ipHE_LIKE ) 00127 { 00128 fprintf(ioQQQ," Sorry, the DAMPING option is not implemented for the he-like sequence.\n"); 00129 cdEXIT(EXIT_FAILURE); 00130 } 00131 00132 /* turn off absorption due to Ly alpha damping wings */ 00133 hydro.DampOnFac = 0.; 00134 } 00135 00136 else if( nMatch("DIEL",chCard) ) 00137 { 00138 if( ipISO == ipH_LIKE ) 00139 { 00140 fprintf(ioQQQ," Sorry, but dielectronic recombination onto the h-like sequence is not possible.\n"); 00141 cdEXIT(EXIT_FAILURE); 00142 } 00143 00144 /* This sets which set of data to use for dielectronic recombination. */ 00145 if( nMatch(" OFF",chCard) ) 00146 { 00147 iso.lgDielRecom[ipISO] = false; 00148 } 00149 else 00150 iso.lgDielRecom[ipISO] = true; 00151 } 00152 00153 else if( nMatch("LEVE",chCard) ) 00154 { 00155 /* the number of levels read in is n, the principal quantum number 00156 * only lines with upper levels less than n will be printed */ 00157 00158 /* number of levels for iso-sequence */ 00159 /* there are two options here, 00160 * when keyword ELEMENT appears, scan off element name and change levels only for 00161 * that one. 00162 * when there is no ELEMENT then set all in iso to same number */ 00163 00164 /* lgHydroMalloc is false at start of calculation, set true when space 00165 * allocated for the hydrogen and helium lines. Once done we must ignore all 00166 * future changes in the number of levels */ 00167 if( nMatch("PRIN",chCard) ) 00168 { 00169 /* only print - do not change levels */ 00170 iso.lgPrintNumberOfLevels = true; 00171 } 00172 else if( !lgHydroMalloc ) 00173 { 00174 i = 5; 00175 numLevels = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00176 00177 if( lgEOL ) 00178 { 00179 /* must be a number, or a key, either large (or limit) or small 00180 * if no number but LARGER or SMALL then set to preset number */ 00181 if( nMatch("LARG",chCard) ) 00182 { 00183 /* there is no real limit, but this includes all levels with tabulated rec coefficient */ 00184 numLevels = RREC_MAXN; 00185 } 00186 00187 /* this is small or compact keyword */ 00188 else if( nMatch("SMAL",chCard) || nMatch("COMP",chCard) ) 00189 { 00190 if( ipISO == ipH_LIKE ) 00191 numLevels = 5; 00192 else if( ipISO == ipHE_LIKE ) 00193 numLevels = 3; 00194 } 00195 00196 else 00197 /* punch out if no number */ 00198 NoNumb(chCard); 00199 } 00200 else if( ( numLevels < 3 ) && !nMatch("COLL",chCard) && ipISO == ipH_LIKE ) 00201 { 00202 /* only enforce this limit if not working on collapsed levels */ 00203 fprintf( ioQQQ, " cannot have fewer than 3 levels, the requested number was %li\n" , 00204 numLevels ); 00205 fprintf( ioQQQ, " Sorry.\n" ); 00206 cdEXIT(EXIT_FAILURE); 00207 } 00208 else if( ( numLevels < 3 ) && !nMatch("COLL",chCard) && ipISO == ipHE_LIKE) 00209 { 00210 /* only enforce this limit if not working on collapsed levels */ 00211 fprintf( ioQQQ, " cannot have fewer than 3 levels, the requested number was %li\n" , 00212 numLevels ); 00213 fprintf( ioQQQ, " Sorry.\n" ); 00214 cdEXIT(EXIT_FAILURE); 00215 } 00216 00217 if( ipISO == ipH_LIKE && numLevels > NHYDRO_MAX_LEVEL-2 ) 00218 { 00219 fprintf( ioQQQ, " Not possible to set nhlvl to >NHYDRO_MAX_LEVEL-2= %i\n", 00220 NHYDRO_MAX_LEVEL-2 ); 00221 fprintf( ioQQQ, " change NHYDRO_MAX_LEVEL\n"); 00222 cdEXIT(EXIT_FAILURE); 00223 } 00224 00225 /* this is check that alpha transition of highest level is within energy 00226 * bounds of continuum */ 00227 if( ipISO == ipH_LIKE && 00228 ( 2. / POW3((double)numLevels) < rfield.emm ) ) 00229 { 00230 fprintf( ioQQQ, " Not possible to set iso.numLevels_max[ipH_LIKE][ipHYDROGEN] to such a high value, since " 00231 "alpha transition not within energy bounds of code\n"); 00232 00233 fprintf( ioQQQ, " lowest energy is %e and corresponding highest level is %li\n" , 00234 rfield.emm, (long)pow(2./rfield.emm, 0.3333) ); 00235 cdEXIT(EXIT_FAILURE); 00236 } 00237 00238 if( nMatch("COLL",chCard) ) 00239 { 00240 if( numLevels < 1 ) 00241 { 00242 fprintf( ioQQQ, "There must be at least one collapsed level.\n"); 00243 cdEXIT(EXIT_FAILURE); 00244 } 00245 00246 /* set collapsed levels for entire sequence or just one element */ 00247 if( nMatch("ELEM",chCard) ) 00248 { 00249 /* check which element is on the line, nelem is atomic number on C scale */ 00250 nelem = GetElem(chCard); 00251 iso.nCollapsed_max[ipISO][nelem] = numLevels; 00252 iso_update_num_levels( ipISO, nelem ); 00253 } 00254 else 00255 { 00256 /* keyword element did not appear, so set all iso species to this number */ 00257 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00258 { 00259 iso.nCollapsed_max[ipISO][nelem] = numLevels; 00260 iso_update_num_levels( ipISO, nelem ); 00261 } 00262 } 00263 } 00264 00265 /* we now have the desired number of levels, and it has been fully qualified. 00266 * now check whether ELEMENT keyword is on line, if not then set all to 00267 * the number, if so then only element in use */ 00268 else if( nMatch("ELEM",chCard) ) 00269 { 00270 /* check which element is on the line, nelem is atomic number on C scale */ 00271 nelem = GetElem(chCard); 00272 iso.n_HighestResolved_max[ipISO][nelem] = numLevels; 00273 iso_update_num_levels( ipISO, nelem ); 00274 } 00275 else 00276 { 00277 /* keyword element did not appear, so set all species to this number */ 00278 for( nelem=ipISO; nelem<LIMELM; ++nelem ) 00279 { 00280 iso.n_HighestResolved_max[ipISO][nelem] = numLevels; 00281 iso_update_num_levels( ipISO, nelem ); 00282 } 00283 } 00284 } 00285 } 00286 00287 else if( nMatch("ERRO",chCard) && nMatch("GENE" , chCard ) ) 00288 { 00289 /* Rates will be modified by a randomly generated error that falls within 00290 * the range specifically set for each rate (or set of rates). */ 00291 iso.lgRandErrGen[ipISO] = true; 00292 i = 5; 00293 iso.modelRank[ipISO] = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00294 00295 iso.modelRank[ipISO] = MAX2( 0, iso.modelRank[ipISO] ); 00296 if( lgEOL ) 00297 /* Changed to avoid lint complaint */ 00298 /* iso.modelRank[ipISO] = (unsigned)time(NULL); */ 00299 iso.modelRank[ipISO] = abs((int)time(NULL)); 00300 00301 /* this allows a seed that is dependent upon the processor rank 00302 * in a parallel run. */ 00303 /* We add 2 so that the seed is always greater than 1, which would reset the generator. */ 00304 init_genrand( (unsigned)iso.modelRank[ipISO] + 2); 00305 } 00306 00307 else if( nMatch(" FSM",chCard) ) 00308 { 00309 if( ipISO == ipH_LIKE ) 00310 { 00311 fprintf(ioQQQ," Sorry, but fine-structure mixing can only be implemented for the He-like sequence.\n"); 00312 cdEXIT(EXIT_FAILURE); 00313 } 00314 00315 /* turn on fine structure mixing of spontaneous decays. 00316 * >>refer Helike FSM Bauman, R., MacAdams, K., and Ferland, G. (2003). */ 00317 if( nMatch(" OFF",chCard) ) 00318 iso.lgFSM[ipISO] = false; 00319 else 00320 iso.lgFSM[ipISO] = true; 00321 } 00322 00323 else if( nMatch("GBAR",chCard) ) 00324 { 00325 if( ipISO == ipH_LIKE ) 00326 { 00327 fprintf(ioQQQ," Sorry, the GBAR option is only implemented for the He-like sequence.\n"); 00328 cdEXIT(EXIT_FAILURE); 00329 } 00330 00331 /* the HEGBAR command - to change cs of higher levels */ 00332 /* first turn all off, one will be turned back on */ 00333 iso.lgCS_Vriens[ipISO] = false; 00334 iso.lgCS_None[ipISO] = false; 00335 iso.nCS_new[ipISO] = false; 00336 00337 /* now turn one on */ 00338 if( nMatch("VRIE",chCard) ) 00339 { 00340 /* option to change how collisions for unknown levels affect things */ 00341 iso.lgCS_Vriens[ipISO] = true; 00342 } 00343 else if( nMatch(" NEW",chCard) ) 00344 { 00345 /* option to change how collisions for unknown levels affect things */ 00346 i = 5; 00347 iso.nCS_new[ipISO] = (int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00348 00349 /* there are two options, 1 and 2, for which fit - default (no number) 00350 * will be 1, the broken power law fit */ 00351 if( lgEOL ) 00352 iso.nCS_new[ipISO] = 1; 00353 00354 ASSERT( iso.nCS_new[ipISO] ); 00355 } 00356 else if( nMatch(" OFF",chCard) ) 00357 { 00358 /* option to change how collisions for unknown levels affect things */ 00359 iso.lgCS_None[ipISO] = true; 00360 } 00361 else 00362 { 00363 fprintf( ioQQQ, " needs parameter\n" ); 00364 cdEXIT(EXIT_FAILURE); 00365 } 00366 } 00367 00368 00369 else if( nMatch("LYMA",chCard) ) 00370 { 00371 if( ipISO == ipH_LIKE && nMatch("PUMP",chCard) ) 00372 { 00373 /* >>chng 05 jul 08, separate out Lyman pump commands */ 00374 if( nMatch(" OFF",chCard) ) 00375 { 00376 /* option to turn off all continuum pumping of Lyman lines */ 00377 hydro.lgLymanPumping = false; 00378 } 00379 else if( nMatch("SCALE",chCard) ) 00380 { 00381 /* multiplicative factor for all continuum pumping of H I Lyman lines, 00382 * account for possible emission in the line - only affects H I 00383 * not entire H-like iso sequence */ 00384 i = 5; 00385 hydro.xLymanPumpingScaleFactor = 00386 (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00387 /* scale factor is log if negative */ 00388 if( hydro.xLymanPumpingScaleFactor < 0. ) 00389 hydro.xLymanPumpingScaleFactor = 00390 (realnum)pow((realnum)10.f , hydro.xLymanPumpingScaleFactor ); 00391 } 00392 else 00393 { 00394 fprintf(ioQQQ," Sorry, I didn\'t recognize an option on this ATOM H-LIKE LYMAN PUMP command.\n"); 00395 fprintf(ioQQQ," The options are \" OFF\", and \"SCALE\".\n"); 00396 cdEXIT(EXIT_FAILURE); 00397 } 00398 } 00399 else 00400 { 00401 /* option to set number of "extra" Lyman lines, used for optical depths only */ 00402 i = 5; 00403 iso.nLyman[ipISO] = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL); 00404 if( lgEOL ) 00405 NoNumb(chCard); 00406 } 00407 } 00408 00409 /* don't interpolate on look-up tables but compute recombination on the fly instead */ 00410 else if( nMatch("RECO",chCard) && 00411 nMatch(" NO ",chCard) && nMatch("INTE",chCard) ) 00412 { 00413 /* flag set by atom xx-like no recombination interp command, 00414 * says to generate recombination coefficients 00415 * on the fly */ 00416 iso.lgNoRecombInterp[ipISO] = true; 00417 } 00418 00419 else if( nMatch("REDI",chCard) ) 00420 { 00421 int ipRedis=0; 00422 /* there are three functions, PRD_, CRD_, and CRDW, 00423 * representing partial redistribution, 00424 * complete redistribution (doppler core only, no wings) 00425 * and complete with wings */ 00426 /* partial redistribution */ 00427 if( nMatch(" PRD",chCard) ) 00428 { 00429 ipRedis = ipPRD; 00430 } 00431 /* complete redistribution no wings */ 00432 else if( nMatch(" CRD",chCard) ) 00433 { 00434 ipRedis = ipCRD; 00435 } 00436 /* complete redistribution with wings */ 00437 else if( nMatch("CRDW",chCard) ) 00438 { 00439 ipRedis = ipCRDW; 00440 } 00441 00442 /* if not SHOW option (handled below) then we have a problem */ 00443 else if( !nMatch("SHOW",chCard) ) 00444 { 00445 fprintf(ioQQQ," There should have been a second keyword on this command.\n"); 00446 fprintf(ioQQQ," Options are _PRD, _CRD, CRDW (_ is space). Sorry.\n"); 00447 cdEXIT(EXIT_FAILURE); 00448 } 00449 00450 /* resonance lines - not Lya*/ 00451 if( nMatch("ALPH",chCard) ) 00452 { 00453 iso.ipLyaRedist[ipISO] = ipRedis; 00454 } 00455 /* Lya itself */ 00456 if( nMatch("RESO",chCard) ) 00457 { 00458 iso.ipResoRedist[ipISO] = ipRedis; 00459 } 00460 /* subordinate lines */ 00461 else if( nMatch("SUBO",chCard) ) 00462 { 00463 iso.ipSubRedist[ipISO] = ipRedis; 00464 } 00465 /* the show option, say what we are assuming */ 00466 else if( nMatch("SHOW",chCard) ) 00467 { 00468 fprintf(ioQQQ," Ly a is "); 00469 if( iso.ipLyaRedist[ipISO] ==ipCRDW ) 00470 { 00471 fprintf(ioQQQ,"complete redistribution with wings\n"); 00472 } 00473 else if( iso.ipLyaRedist[ipISO] ==ipCRD ) 00474 { 00475 fprintf(ioQQQ,"complete redistribution with core only.\n"); 00476 } 00477 else if( iso.ipLyaRedist[ipISO] ==ipPRD ) 00478 { 00479 fprintf(ioQQQ,"partial redistribution.\n"); 00480 } 00481 else if( iso.ipLyaRedist[ipISO] ==ipLY_A ) 00482 { 00483 fprintf(ioQQQ,"special Lya.\n"); 00484 } 00485 else 00486 { 00487 fprintf(ioQQQ," PROBLEM Impossible value for iso.ipLyaRedist.\n"); 00488 TotalInsanity(); 00489 } 00490 00491 fprintf(ioQQQ," Other %s resonance lines are ", 00492 elementnames.chElementSym[ipISO] ); 00493 00494 if( iso.ipResoRedist[ipISO] ==ipCRDW ) 00495 { 00496 fprintf(ioQQQ,"complete redistribution with wings\n"); 00497 } 00498 else if( iso.ipResoRedist[ipISO] ==ipCRD ) 00499 { 00500 fprintf(ioQQQ,"complete redistribution with core only.\n"); 00501 } 00502 else if( iso.ipResoRedist[ipISO] ==ipPRD ) 00503 { 00504 fprintf(ioQQQ,"partial redistribution.\n"); 00505 } 00506 else 00507 { 00508 fprintf(ioQQQ," PROBLEM Impossible value for iso.ipResoRedist.\n"); 00509 TotalInsanity(); 00510 } 00511 00512 fprintf(ioQQQ," %s subordinate lines are ", 00513 elementnames.chElementSym[ipISO] ); 00514 00515 if( iso.ipSubRedist[ipISO] ==ipCRDW ) 00516 { 00517 fprintf(ioQQQ,"complete redistribution with wings\n"); 00518 } 00519 else if( iso.ipSubRedist[ipISO] ==ipCRD ) 00520 { 00521 fprintf(ioQQQ,"complete redistribution with core only.\n"); 00522 } 00523 else if( iso.ipSubRedist[ipISO] ==ipPRD ) 00524 { 00525 fprintf(ioQQQ,"partial redistribution.\n"); 00526 } 00527 else 00528 { 00529 fprintf(ioQQQ," PROBLEM Impossible value for iso.ipSubRedist.\n"); 00530 TotalInsanity(); 00531 } 00532 } 00533 else 00534 { 00535 fprintf(ioQQQ," here should have been another keyword on this command.\n"); 00536 fprintf(ioQQQ," Options are ALPHA, RESONANCE, SUBORDINATE. Sorry.\n"); 00537 cdEXIT(EXIT_FAILURE); 00538 } 00539 } 00540 00541 /* which type of matrix solution, populations or lowt */ 00542 else if( nMatch("MATR",chCard) ) 00543 { 00544 if( nMatch(" POP",chCard) ) 00545 { 00546 strcpy( iso.chTypeAtomSet[ipISO] , "POPU" ); 00547 } 00548 else if( nMatch(" LOW",chCard) ) 00549 { 00550 strcpy( iso.chTypeAtomSet[ipISO] , "LOWT" ); 00551 } 00552 else 00553 { 00554 fprintf( ioQQQ, " There should have been a keyword on the MATRIX option.\n"); 00555 fprintf( ioQQQ, " The keywords are POPulations, and LOWte.\n"); 00556 cdEXIT(EXIT_FAILURE); 00557 } 00558 } 00559 00560 else if( nMatch("TOPO",chCard) ) 00561 { 00562 if( ipISO == ipH_LIKE ) 00563 { 00564 fprintf( ioQQQ, " The behavior of the TOPOFF option for the H-like sequence has changed since version 07.02.\n"); 00565 fprintf( ioQQQ, " Originally, it accepted a number and distributed \"topoff\" equally into that many levels at the top of the model.\n"); 00566 fprintf( ioQQQ, " Now, the option works the same as in the He-like sequence, which is to turn off topoff altogether.\n" ); 00567 } 00568 00569 if( nMatch(" OFF",chCard) ) 00570 { 00571 iso.lgTopoff[ipISO] = false; 00572 fprintf( ioQQQ, "ISO %li TOPOFF is OFF\n", ipISO ); 00573 } 00574 else 00575 iso.lgTopoff[ipISO] = true; 00576 } 00577 00578 00579 else 00580 { 00581 fprintf( ioQQQ, " There should have been a keyword - STOP\n" ); 00582 cdEXIT(EXIT_FAILURE); 00583 } 00584 return; 00585 }